Algorithmen Teil VIII: Ein Bug kommt selten allein...

WindowTitle "Datenstruktur DEQue (DoubleEnded Queue): Die beidseitige Warteschlange (aka Doppelstack)"
'(CL) CopyLeft 2014-08 by P. Specht, Vienna/Austria 
' Keine wie auch immer geartete Gewähr, Verwendung auf alleinige Gefahr der AnwenderInnen! 
WindowStyle 24:CLS 'Window 0,0-%maxx,%maxy-40
Font 2
Print "\n============================================================"
print "                    KLEINES DEQue-TUTORIAL                  "
print "  PROCs are Freeware,  Tutorial=CopyLeft by P. Specht       "
Print "============================================================"
Print "\n Beim Umsetzen z.B. des Melkman Convex Hull-Algorithmus [Melkman 1987]    "
Print " begegnet man einer Datenstruktur namens DEQue (spr.'Deck' Double Ended Queue)"
print " Während 'Queue' im Englischen 'Warteschlange' bedeutet (aka. FIFO-Speicher)"
print " handelt es sich beim DEQue um eine Warteschlange, bei der zB. Privatpatien-"
print " ten vorgezogen werden können, aber auch Leute die frustrierend lange warten"
print " sich hinten wieder entfernen (oder die Schlange wechseln) können.          "
print "\n Auf Datenhaltungsstrukturen bezogen bedeutet die Realisierung, daß entweder"
print " Arrayelemente physisch stellenweise vor- bzw. zurückverschoben werden (hier"
print " im Programm realisiert, relativ langsam), oder vor-&rück-verlinkte Listen  "
print " die entsprechende Funktionalität bereitstellen (Aufwendig, am besten als   "
print " Memory-Verwaltungssystem auf maschinennaher Ebene, Betriebssystemabhängig)."
print "\n Hinweis: Einfache Warteschlangen benötigten nur die Hälfte der PROCs hier! "
print "\n Der Begriff TOP bedeutet hier 'OBEN', der Begriff Btm=Bottom ='UNTEN'      "
print "\n Quelle: http://de.wikipedia.org/wiki/Deque                          [Taste]"
Font 0
waitinput
Cls


Declare Wert&,DEQ&[],numw&
numw&=4 'temporary numwidth, mit der das DEQue zu Testzwecken ausgegeben wird


print "\n   MakeConstDEQue DEQ&[],11,-7" 
MakeConstDEQue DEQ&[],11,-7
ShowDEQue DEQ&[],numw&
Waitinput
print "\n   MakeConstDEQue DEQ&[],11,88"
MakeConstDEQue DEQ&[],11,88
ShowDEQue DEQ&[],numw&
Waitinput
print "\n   MakeConstDEQue DEQ&[],11     <<< Ohne Wertangabe wird NICHT überschrieben!"
MakeConstDEQue DEQ&[],11
ShowDEQue DEQ&[],numw&
Waitinput
print "\n   MakeIndexDEQue DEQ&[],7"
MakeIndexDEQue DEQ&[],7
ShowDEQue DEQ&[],numw&
Waitinput
print "\n   MakeRevertIndexDEQue DEQ&[],11"
MakeRevertIndexDEQue DEQ&[],11
ShowDEQue DEQ&[],numw&
Waitinput


skip:
CLS
print "\n   Clear DEQ&[]"
Clear DEQ&[]
ShowDEQue DEQ&[],numw&
waitinput
Print "\n   Mehrfaches   PutTop(Wert&,DEQ&[])"
DEQ&[]=PutTop(555,DEQ&[])
DEQ&[]=PutTop(111,DEQ&[])
DEQ&[]=PutTop(333,DEQ&[])
DEQ&[]=PutTop(777,DEQ&[])
ShowDEQue DEQ&[],numw&
Waitinput
Print "\n   Mehrfaches PopTop(DEQ&[])   <<< mehr als Werte in der Schlange stehen)"
font 2:print " Top=";PopTop(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
Waitinput
font 2:print " Top=";PopTop(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
Waitinput
font 2:print " Top=";PopTop(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
Waitinput
font 2:print " Top=";PopTop(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
print "   Magic Error-Wert =";int(999999999);" (= minus 9 mal die 9)"
font 2:print " Top=";PopTop(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
font 2:print " Top=";PopTop(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
CLS
Print "\n   Clear DEQ&[]"
Clear DEQ&[]
ShowDEQue DEQ&[],numw&
print
waitinput
print "\n  Mehrfaches  PutBtm(Wert&,DEQ&[])" 
DEQ&[]=PutBtm(666,DEQ&[])
ShowDEQue DEQ&[],numw&
waitinput
DEQ&[]=PutBtm(444,DEQ&[])
ShowDEQue DEQ&[],numw&
waitinput
DEQ&[]=PutBtm(222,DEQ&[])
ShowDEQue DEQ&[],numw&
waitinput


print "\n   Mehrfaches  Wert&=PopBtm(DEQ&[])"
font 2:print " Btm=";PopBtm(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
font 2:print " Btm=";PopBtm(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
font 2:print " Btm=";PopBtm(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
font 2:print " Btm=";PopBtm(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
font 2:print " Btm=";PopBtm(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
CLS
Print "\n   Clear DEQ&[]"
CLEAR DEQ&[]
ShowDEQue DEQ&[],numw&
waitinput
Print "\n   Mehrfaches  PutTop(.) / PutBtm(),  wild durcheinander"
DEQ&[]=PutTop(555,DEQ&[])
DEQ&[]=PutBtm(444,DEQ&[])
DEQ&[]=PutTop(666,DEQ&[])
ShowDEQue DEQ&[],numw&
waitinput
Print "\n   Mehrfaches  PopTop / PopBtm,  wild durcheinander"
font 2:print " Top=";PopTop(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
font 2:print " Btm=";PopBtm(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
font 2:print " Btm=";PopBtm(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
font 2:print " Top=";PopTop(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
CLS
Print "\n   Clear DEQ&[]"
CLEAR DEQ&[]
ShowDEQue DEQ&[],numw&
waitinput
Print "\n   Mehrfaches  PutTop(.) / PutBtm(),  wild durcheinander"
DEQ&[]=PutTop(555,DEQ&[])
DEQ&[]=PutBtm(444,DEQ&[])
DEQ&[]=PutTop(666,DEQ&[])
ShowDEQue DEQ&[],numw&
waitinput
font 2:print " CopyTop=";CopyTop(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
waitinput
font 2:print " CopyBtm=";CopyBtm(DEQ&[]):font 0
ShowDEQue DEQ&[],numw&
Print "\n\n============================================================"
print "             ENDE DES KLEINEN DEQue-TUTORIALS               "
print "                                          Gruss, P. Specht  "
Print "============================================================"
waitinput
End






proc MakeIndexDEQue :parameters q&[],n&
  SetSize q&[],n&
  q&[]=&index+1
return q&[]
endproc


proc MakeRevertIndexDEQue :parameters q&[],n&
SetSize q&[],n&
q&[]=n&-&index
return q&[]
endproc


proc MakeConstDEQue :parameters q&[],n&,wert&
SetSize q&[],n&
q&[]=Wert&
return q&[]
endproc


Proc PopBtm :parameters q&[]
'pops bottom and returns popped bottom element
  var sz&=sizeof(q&[])
  if sz&>0
    var wert&=q&[0]
    q&[]=q&[&index+(&index<(sz&-1))]
    setsize q&[],sizeof(q&[])-(sz&>0)
    return wert&
  else
    return int(-999999999)
  endif
endproc


proc PopTop :parameters q&[]
declare wert&
' pops and returns the popped top element
var sz&=sizeof(q&[]) 
if sz&>0
   wert&=q&[sz&-1]
else
   wert&=int(-999999999)
endif
setsize q&[],sz&-(sz&>0)
return wert&
endproc


Proc PutTop :parameters Wert&,q&[]
  var sz&=sizeof(q&[])
  inc sz&
  setsize q&[],sz&
  q&[sz&-1]=Wert&
endproc


proc PutBtm :parameters Wert&,q&[]
declare tmp&[]
  var sz&=sizeof(q&[])
  inc sz&
  setsize q&[],sz&
  tmp&[]=q&[]
  q&[]=tmp&[&index-(&index>0)]
  clear tmp&[]
  q&[0]=Wert&
endproc


Proc CopyBtm :parameters q&[]
' Returns copy of bottom element
  var sz&=sizeof(q&[])
  if sz&>0
    return q&[0]
  else
    return int(-999999999)
  endif
endproc


proc CopyTop :parameters q&[]
declare wert&
' returns copy of top element
var sz&=sizeof(q&[]) 
if sz&>0
   wert&=q&[sz&-1]
else
   wert&=int(-999999999)
endif
return wert&
endproc


proc ShowDEQue :parameters q&[],nw&
var tmp&=get("numwidth"):set("numwidth",nw&)
  var sz&=sizeof(q&[])
  print tab(2);"n=";sz&;":  ";
  whileloop 0,sz&-1
    print q&[&Loop],  
  endwhile
  print
set("numwidth",tmp&)
endproc

proc CopyFromDEQue :parameters q&[],n&
declare wert&
' returns copy of top element
var sz&=sizeof(q&[]) 
if (sz&>0) and (n&>=0) and (n&<sz&)
   wert&=q&[n&]
else
   wert&=int(-999999999)
endif
return wert&
endproc


Proc CmpInDEQue :parameters q&[],n&,m&
'Compares two elements in DEQue
var sz&=sizeof(q&[]) 
if (sz&>0) and (n&>=0) and (n&<sz&) and (m&>=0) and (m&<sz&)
   return int((q&[n&]>q&[m&])-(q&[n&]<q&[m&]))
else
   return int(-999999999)
endif
Endproc


proc SwapInDEQue :parameters q&[],n&,m&
'Swaps two elements in Queue
var sz&=sizeof(q&[])
if (sz&>0) and (n&>=0) and (n&<sz&) and (m&>=0) and (m&<sz&)
   var tmp&=q&[n&]:q&[n&]=q&[m&]:q&[m&]=tmp&
   return int(n&<>m&)
else
   return int(-999999999)
endif
endproc


proc DelFromDEQue :parameters q&[],n&
'Delets member n from DEQue
var sz&=sizeof(q&[]) 
declare wert&
if (sz&>0) and (n&>0) and (n&<sz&) 
   wert&=q&[n&]
   q&[]=q&[&index+(&index>=n&)*(&index<(sz&-1))]
   dec sz&
   setsize q&[],sz&
else
   wert&=int(-999999999)
endif
return wert&
endproc

Print "  Wert& = CopyFromDEQue(DEQ&[],1) "
font 2:print " n.Element=";CopyInDEQue(DEQ&[],1):font 0
ShowDEQue DEQ&[],numw&
waitinput


Print "  Wert& = CopyFromDEQue(DEQ&[],11111) "
font 2:print " n.Element=";CopyInDEQue(DEQ&[],11111):font 0
ShowDEQue DEQ&[],numw&
waitinput


CLS


print "\n (n>m)? = CmpInDEQue(q&[],0,1)"
font 2:print " LT-1 EQ0 GT1 =";CmpInDEQue(DEQ&[],0,1):font 0
ShowDEQue DEQ&[],numw&
waitinput


Print "\n  (n=m)? = SwapInDEQue(q&[],0,1)"
font 2:print " (N=M)=";SwapInDEQue(DEQ&[],0,1):font 0
ShowDEQue DEQ&[],numw&
waitinput


print "\n (n>m)? = CmpInDEQue(q&[],0,1)"
font 2:print " LT-1 EQ0 GT1 =";CmpInDEQue(DEQ&[],0,1):font 0
ShowDEQue DEQ&[],numw&
waitinput


Print "    ok = DelFromDEQue(DEQ&[],1)"
DEQ&[]=DelFromDEQue(DEQ&[],1)
ShowDEQue DEQ&[],numw&
waitinput

WindowTitle upper$(" P o s i t i o n   d e r   P l a n e t e n   b e s t i m m e n ")+" - Achtung, Daten aus 1980!"
' Source: http://jean-pierre.moreau.pagesperso-orange.fr/Basic/planets_bas.txt
' Translated to XProfan-11 (CL) CopyLeft 2014-08 by P. Specht, Vienna
' No warranty whatsoever! Keine wie auch immer geartete Gewähr!
' Bahndaten aus der Zeit vor 1982, zu verifizieren!
'
'*****************************************************
'*              Position der Planeten                *
'* ------------------------------------------------- *
'* This program calculates the ephemeris of a planet *
'* in our solar system to adjust an equatorial tele- *
'* scope (Pgm für aequatorial-montierte Teleskope    *
'* ------------------------------------------------- *
'* Ref.: "Mathematiques par l'informatique indivi-   *
'*        duelle - Programmes en BASIC, MASSON,      *
'*        Paris, 1982, p 100 - 105" .                *
'* ------------------------------------------------- *
'* Inputs:                                           *
'*   Date: day,month,year                            *
'*   Hour UT: hour                                   *
'*   Planet number: 1 to 8                           *
'*   (Mercury: 1 Venus : 2 Mars   : 4 Jupiter: 5     *
'*    Saturn : 6 Uranus: 7 Neptune: 8)               *
'* Outputs:                                          *
'*   Ascent in hours (0 to 24)                       *
'*   Declination in deg. (-90 to 90 North or South)  *
'*                                                   *
'* SAMPLE RUN:                                       *
'* (Find ascent and declination of planet Mars on    *
'*  March 10th, 1982 at 6h UT)                       *
'                                                    *
' Hinweis von P.Specht:
' Am 10. März 1982 standen Merkur, Venus, Erde, 
' Mars, Jupiter, Saturn und Pluto alle in einer Linie
' und auf der selben Seite der Sonne innerhalb eines 
' Winkel von 95 Grad. Die prophezeiten Katastrophen,
' etwa Erdbeben in Kalifornien sind gsd ausgeblieben.
' http://www.borderlands.de/net_pdf/NET1107S30-34.pdf
'
'* Date (D,M,Y): 10,3,1982                           *
'* Hour UT: 6                                        *
'* Mercury: 1 Venus : 2 Mars   : 4 Jupiter: 5        *
'* Saturn : 6 Uranus: 7 Neptune: 8                   *
'* Planet number: 4                                  *
'*                                                   *
'* Ascent     : 13 H  8 MN                           *
'* Declination:  3 ø 45 MN S                         *
'*                                                   *
'* This program allows building the following table: *
'*                                                   *
'* Date: March 10th, 1982 at 6H UT.                  *
'*                                                   *
'*       Planet     Ascent     Declination           *
'*       Mercury    21 H 51    14 ø 45 mn S          *
'*       Venus      20 H 26    14 ø 57 mn S          *
'*       Mars       13 H 08     3 ø 45 mn S          *
'*       Jupiter    14 H 32    13 ø 30 mn S          *
'*       Saturn     13 H 22     5 ø 42 mn S          *
'*                                                   *
'*             English Version By J-P Moreau, Paris. *
'*                        (www.jpmoreau.fr)          *
'*****************************************************
WindowStyle 24:Window 0,0-%maxx,%maxy:font 2
declare i&,j&,jj&,n&,m&,mm&,aa!,h!,t!,xl!,o!,p!,e!,xi!,xa!,xm!
declare u!,r!,v!,q!,xj!,x!,y!,z!,G!,pi!,xx!,yy!,zz!,d!,b!,c!
declare A![9,8],B![12],A$,data$
'Init planetary constants (fill table A by columns)
data$="4.01166, 0.071425454, 1.32493, 0.000000742289, 0.823045, "+\
"0.000000566185, 0.205615, 0.122225, 0.387099, "+\
"3.60861, 0.027963119, 2.271616, 0.00000065572, 1.32291,"+\
"0.000000436681, 0.006816, 0.0592301, 0.723332,"+\
"1.72727, 0.0172028, 1.76688, 0.000000818559, 0,"+\
"0, 0.016751, 0, 1,"+\
"2.17756, 0.0091467658, 5.83378, 0.000000879297, 0.851616, "+\
"0.000000371232, 0.093309, 0.0322939, 1.5236,"+\
"4.68279, 0.00145099, 0.2289, 0.000000857, 1.73578,"+\
"0.000000482933, 0.048376, 0.0228418, 5.202799,"+\
"4.8567, 0.00058484, 1.5974, 0.000000412, 1.96856,"+\
"0.000000417308, 0.054311, 0.0435026, 9.552098, "+\
"4.3224, 0.000205424, 2.9523, 0.000000762, 1.2825, "+\
"0.000000238237, 0.047319, 0.013482, 19.21694, "+\
"1.5223, 0.000105061, 0.7637, 0.000000393, 2.28102,"+\
"0.00000052517, 0.008262, 0.0310536, 30.12912"
'=Spaltenweise! Schnell (=in der inneren Schleife) läuft daher der ZEILEN-Index i&!)
whileloop 8:j&=&Loop:whileLoop 9:i&=&loop
  A![i&,j&]=val(substr$(data$,(i&+(j&-1)*9),","))
' print i&;",";j&;"=",int(1+(i&-1)+(j&-1)*9);":",format$("%g",a![i&,j&]),"  ",:waitinput
endwhile :endwhile ' Basis [1,1] wird hier beibehalten
'Calendar constants (Ermittlung der Tagesnummer im Jahr)
'Sinnstörender Komma-Fehler im Originaltext!
'Schaltjahre sind hier aufgeteilt: Die Erde selbst schaltet ja nicht ;-)
whileloop 12 :B![&Loop]=\
val(substr$("0,31,59.25,90.25,120.25,151.25,181.25,212.25,243.25,273.25,304.25,334.25",&Loop,","))
endwhile
print "\n\n\n "
'print " Please input date (DD, MM, YYYY) and hour (GMT) of interest:"
print " Bitte das interessierende Datum (DD, MM, YYYY) und die UTC-Stunde (GMT) eingeben:"
print "\n  Tag xx:",:input j&
print "  Monat xx:",:input mm&
print " Jahr xxxx:",:input aa!
print "\n Stunde UT = GMT(Greenwich): ",:input h!


'Calculate time t
t!=365.25*(aa!-1901)+B![mm&]+j& 'Tagesnummer seit 1901
t!=INT(t!)+h!/24                'Stunden in Dezimalform dazu
'calculate earth coordinates
j&=3 : gosub "S1000"    'planet #3 coordinates
g!=x! : h!=y!           'save earth coordinates


Nochmal:
print
print " Merkur/y=1, Venus=2,((Earth=3)), Mars=4, Jupiter=5, Saturn=6, Uranus=7, Neptun=8"
print " Please Input Planet Number: ", :input n&
print 
'calculate coordinates of planet #n
j&=n& : gosub "S1000"
'calculate geocentric equatorial coordinates
x!=x!-g! : y!=y!-h!
t!=y!*0.917484-z!*0.397772
z!=y!*0.397772+z!*0.917484 : y!=t!
'calculate ascent and declination
pi!=pi() '=4*arctan(1)
xx!=x! : yy!=y! : gosub "S2000"  'aa=ATN2(y,x) ' Simuliert echte Atn2( , )-Funktion
aa!=zz!
yy!=z! : xx!=SQRT(x!*x!+y!*y!) : gosub "S2000" ' d=ATN2(z, SQR(x*x+y*y) )
d!=zz!
'conversion
aa!=aa!*12/pi!
case aa!<0 :aa!=24+aa!
h!=INT(aa!) : m&=INT(60*(aa!-h!))
d!=d!*180/pi!
A$="N":case d!<0:A$="S"
d!=ABS(d!) : b!=INT(d!) : c!=INT(60*(d!-b!))
'print results
print
print "\n           Ascension (Aufgangszeit): ";format$("00",h!);" Uhr ";format$("00",m&);" min"
print "\n bei Deklination (Fusspunkt-Winkel): ";
print format$("#00",b!);"° ";format$("00",c!);" ";A$
waitinput
' Auf Grund der unglücklichen Zweitverwendung von Variablen im Original stets Neustart erforderlich!
END


S1000:
'planet coordinates
'calculate planetary constants
  XL!=A![1,j&] + A![2,j&] * t!
  O! =A![3,j&] + A![4,j&] * t!
  P! =A![5,j&] + A![6,j&] * t!
  E! =A![7,j&]
  XI!=A![8,j&]
  XA!=A![9,j&]
'solve Kepler's equation
  xm! = XL!-O! : u!=xm!
  whileloop 100:j&=&Loop
    u! = xm! + E! * sin(u!)
  endwhile


'formula (3) of reference book
 r!=XA!*(cos(u!)-E!)
 v!=XA!*SQRT(1-E!*E!)*sin(u!)
 O!=O!-P! : xm!=sin(O!) : O!=cos(O!)
 q!=sin(P!) : P!=cos(P!)
 xj!=sin(XI!) : XI!=cos(XI!)
 x!=(P!*O!-XI!*q!*xm!)*r!+(-P!*xm!-XI!*q!*O!)*v!
 y!=(q!*O!+XI!*P!*xm!)*r!+(-q!*xm!+XI!*P!*O!)*v!
 z!=xj!*xm!*r!+xj!*O!*v!
return


S2000:
'calculate zz = ATN2(yy,xx)
if ABS(xx!)<val("1e-12")
  if yy!>0 : zz!=Pi!/2
  elseif ABS(yy!)<val("1e-12"): zz!=-PI!/2
  else : zz!=0
  endif
else
  zz!=Arctan(yy!/xx!)
  case xx!<0 :zz!=zz!+Pi!
endif
return
'end of file planets.prf

WindowTitle "Anfangswertproblem einer Differentalgleichung 2. Ordnung lösen nach Størmer-Verlet"
' ---------------------------------------------------------------------------------------------
' Aus Fortran90 nach XProfan-11 übersetzt 2014-08 von P. Specht, Wien (Österreich)
' Keine wie auch immer geartete Gewähr! Verwendung auf alleinige Gefahr de(s|r) Anwender(s|in)!
' Quelle: http://jean-pierre.moreau.pagesperso-orange.fr/Cplus/stormer.txt
'**********************************************************************************************
'*  Differential equation y"=f(x,y,y') by  Størmer's method  *
' (erweiterbar auf allg. Verlet type, z.B. y'''[...]=f(x,y,y',y''[,y''',...])
'* --------------------------------------------------------- *
'* SAMPLE RUN (BEISPIEL)                                     *
'* Löse  y" = 8*y*y / (1+2*x) from x=0 to x=1                *
proc G :parameters x!,y!,z!  '= x, y, y'
  return 8*sqr(y!)/(1+2*x!)
endproc
proc F :parameters x!,y!,z!  '= x, y, y'
  return z!
endproc
'* with initial conditions: x(0)=0, y(0)=1 and y'(0)=-2      *
Declare x![4],y![4],z![4],h!  ' Initial conditions:
x![1]= 0
y![1]= 1
z![1]=-2  ' dy/dt
h!=0.01   ' Step
'*  Compare to exact solution   y = 1/(1+2*x)                *
Proc Fx :parameters x! ' Exact Solution:
  return 1/(1+2*x!)
endproc
'* --------------------------------------------------------- *
'*   Differential equation y"=f(x,y,y') by Stormer's method  *       
'* --------------------------------------------------------- *
'*        X           Y           Y exact         Error      *
'*     0.000       1.000000       1.000000     0.0000000000  *
'*     0.010       0.980392       0.980392     0.0000000001  *
'*     0.020       0.961538       0.961538     0.0000000295  *
'*     0.030       0.943396       0.943396     0.0000000457  *
'*     0.040       0.925926       0.925926     0.0000000974  *
'*     0.050       0.909091       0.909091     0.0000001285  *
'*     ...         ...            ...          ...           *
'*     0.950       0.344866       0.344828     0.0000381695  *
'*     0.960       0.342505       0.342466     0.0000388874  *
'*     0.970       0.340176       0.340136     0.0000396196  *
'*     0.980       0.337878       0.337838     0.0000403406  *
'*     0.990       0.335612       0.335570     0.0000410721  *
'*     1.000       0.333375       0.333333     0.0000418231  *
'*                                                           *
'*   F90 Release By J-P Moreau, Paris (www.jpmoreau.fr).     *
'*************************************************************
WindowStyle 24:font 2:set("decimals",16):Window 0,0-%maxx,%maxy-40
Declare c![4],a1!,a2!,a3!,yex!,k&,er! , x1!,y1!,z1!
a1!= 1.083333333333333333 : a2!=-2*(a1!-1) : a3!= a1!-1
clearclip
print "\n------------------------------------------------------------------------------"
print "       Differentialgleichung y''=f(x,y,y') mit Störmer's Methode lösen "
print "------------------------------------------------------------------------------\n"
print "         X                   Y                  Y exact                Error\n"
putclip " X | Y | Y_exact | Error \n"
yex!=Fx(x![1]):print tab(2);x![1],tab(22);y![1],tab(42);yex!,tab(62);er!
putclip str$(x![1])+"|"+str$(y![1])+"|"+str$(yex!)+"|"+str$(er!)+"\n"
' Runge-Kutta for first 2 steps:
whileloop 1,2:k&=&Loop  
  RK4D2 x![k&],y![k&],z![k&],h!
  x![k&+1]=x1!
  y![k&+1]=y1!
  z![k&+1]=z1!
  yex!=Fx(x![k&+1])
  er!=abs(yex!-y![k&+1])
print tab(2);x![k&+1],tab(22);y![k&+1],tab(42);yex!,tab(62);er!    
putclip str$(x![k&+1])+"|"+str$(y![k&+1])+"|"+str$(yex!)+"|"+str$(er!)+"\n"
endwhile


REPEAT 'Main loop G10:
  whileloop 2,4:k&=&Loop 
    c![k&]=G( x![5-k&],y![5-k&],z![5-k&] )
  endwhile
  y![4]=2*y![3]-y![2]+sqr(h!)*(a1!*c![2]+a2!*c![3]+a3!*c![4])
  x![4]=x![3]+h!
  yex!=Fx(x![4])
  er!=abs(yex!-y![4])
  print tab(2);x![4],tab(22);y![4],tab(42);yex!,tab(62);er!
  putclip str$(x![4])+"|"+str$(y![4])+"|"+str$(yex!)+"|"+str$(er!)+"\n"
  Whileloop 1,3:k&=&Loop 'für nächsten Schritt umkopieren:
    x![k&]=x![k&+1]
    y![k&]=y![k&+1]
  endwhile
  if %csrlin>40:waitinput 4000:cls
    print "\n         X                   Y                  Y-exakt                Error\n"
  endif
UNTIL x![3]>1 'end x value = 1
print "\n EOF"
Print " Resultate in Zwischenablage!"
Waitinput
END


PROC RK4D2 :parameters x!,y!,z!,h! ',x1!,y1!,z1!
declare c1!,d1!,c2!,d2!,c3!,d3!,c4!,d4!
c1!=F(x!,y!,z!)
d1!=G(x!,y!,z!)
c2!=F(x!+h!/2, y!+h!/2*c1!, z!+h!/2*d1!)
d2!=G(x!+h!/2, y!+h!/2*c1!, z!+h!/2*d1!)
c3!=F(x!+h!/2, y!+h!/2*c2!, z!+h!/2*d2!)
d3!=G(x!+h!/2, y!+h!/2*c2!, z!+h!/2*d2!)
c4!=F(x!+h!,   y!+h!*c3!,   z!+h!*d3!  )
d4!=G(x!+h!,   y!+h!*c3!,   z!+h!*d3!  )
x1!=x!+h!
y1!=y!+h!*(c1!+2*c2!+2*c3!+c4!)/6
z1!=z!+h!*(d1!+2*d2!+2*d3!+d4!)/6
endproc
PROGEND
'{=============================================================================================
'    Stormer, aka Newton-Størmer-Verlet leapfrog method
' =============================================================================================
' Namensgebung: Der Algorithmus selbst geht zurück auf Delambre (1791), oftmals wiederentdeckt,
' u.a. 1907 durch Carl Størmer, 1909 durch Cowell und Cromlin, 1967 populär durch Loup Verlet.
' Verbessert die explizite Euler-Methode ohne allzu hohen Mehraufwand. Oft in Physik-Engines
' zu finden. Details siehe http://en.wikipedia.org/wiki/Verlet_integration
' =============================================================================================
' Methode STORMER (aka Newton-Störmer-Verlet-leapfrog method):
' löst Y"=f(x,y,y') mit gegebenen Anfangsbedingungen
' ----------------------------------------------------------------------------------
' Eine Differentialgleichung 2. Ordnung (Y") kann ersetzt werden durch ein
' Differentialgleichungs-SYSTEM bestehend aus 2 Gleichungen 1. Ordnung!
' Wenn also y"=f(x,y,y') gegeben ist mit y(a=Anfang) und y'(a=Anfang), dann kann man
' diese Formel durch Substitution u = y' umformen zu: 
'
'                    | u'=f(x,y,u)
'                    | y'=u
'                    | mit den neuen Anfangsbedingungen y(a), u(a)
'
' Beginnen wir die Betrachtung mit jener speziellen Form der Taylor-Entwicklung, 
' bei der der Rest als Integral dargestellt wird (Wiki: Taylorentwicklung)
'                                       x+h
'           y(x+h) = y(x) + h * y'(x) + INTEGRAL (x+h-t)*y"(t).dt
'                                        x
' Durch Verbringung der Nichtintegral-Terme auf die linke Seite und Aufspaltung
' des 'Restfehler-Integrals' in zwei Teile kann diese Form auch folgender-
' maßen dargestellt werden:
'                                        x+h              x-h
'           y(x+h) - 2y(x) + y(x-h) = INTEGRAL ...  +  INTEGRAL ...
'                                         x                x
' Für das zweite Integral wird die Variable x nun substituiert durch:  u = 2*x-t
' Dann gilt auf Grund der Kettenregel:
'           x-h                              x+h
'         INTEGRAL (x-h-t) * y"(t).dt = - INTEGRAL (x+h-u)*y"(u).du
'            x                                x
' und die Ableitung .du rück-substituiert (zu .-dt) liefert dann:
'                                    x+h
'           du=-dt ==>           = INTEGRAL (x+h-t)*y"(2x-t).dt
'                                     x
' Somit kann man schreiben:              x+h
'           y(x+h) - 2*y(x) + y(x-h) = INTEGRAL (x+h-t)*[y"(t)+y"(2x-t)].dt
'                                         x
' Nun verwenden wir unter Hinweis, daß ja y"(x)=f(x,y,y') gilt, folgenden Ausdruck 
' als Interpolationspolynom:                               
'                             x_n+1    
'     y_n+1 - y_n + y_n-1 = INTEGRAL (x_n+1 - t)*[P(t)+P(2*x_n - 1)].dt  = ...
'                              x_n             
'               x_n+1                                 
'    ...= h^2*INTEGRAL (x_n+1 -t)*[O_00*Div(f_n)+O_11*Div(f_n)+O_22*Div(f_n)].dt ,
'                x_n
'                         1           |(-s) (m)|
' wobei O_mm  = -1^m * INTEGRAL (1-s)*|    +   | .ds
'                         0           |( m) (s)|
'                                               (vgl. auch 'Adams-Bashford-Methode')
' Letzteres bringt uns schließlich zu STORMER's Formeln:
' ==================================================================================
' Explizit formuliert:  y_n+1 -2*y_n + y_n-1 = h^2/12 * [13 * f_n - 2*f_n-1 + f_n-2]
' Das interessierende   y_n+1 = h^2/12 * [13 * f_n - 2*f_n-1 + f_n-2]+ 2*y_n - y_n-1 
' (f_n+1 muss nicht bekannt sein: Prognoserechnung bei aktuell ergänzten Zeitreihen)
' ----------------------------------------------------------------------------------
' oder implizit:        y_n+1 -2*y_n + y_n-1 = h^2/12 * [f_n+1 + 10 * f_n + f_n-1 ]
'}==================================================================================

WindowTitle " Anfangswert-bestimmte Differentialgleichung "+\
"1. Ordnung mit PRÄDIKTOR-KORREKTOR-METHODE lösen"
'{ (T) Translated to XProfan-11 2014-08 by P.Specht, Wien
'***********************************************************
'*        SOLVING DIFFERENTIAL EQUATIONS OF ORDER 1        *
'*             (Prediction-correction method)              *
'* ------------------------------------------------------- *
' Beisp. für eine DGL 1. Ordnung: Etwa das Zerfallsgesetz  *
' d.Radioaktivität. Nachstehend: Willkürliche Testfunktion *
'* ------------------------------------------------------- *
'* Reference:                                              *
'*   "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc  *
'*    DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991"    *
'* F90 version by J-P Moreau, Paris   (www.jpmoreau.fr)    *
'* ------------------------------------------------------- *
'* SAMPLE RUN:                                             *
'* (Integrate y'=(4xy+4x*sqrt(y))/(1+x²) from x=0 to x=10  *
'*  with initial condition y(0)=1)                         *
'*                                                         *
'*  Input x begin: 0                                       *
'*  Input x end  : 10                                      *
'*  Input y begin: 1                                       *
'*  Input number of points: 10                             *
'*  Input finesse: 10                                      *
'*                                                         *
'*         X               Y                               *
'*  -----------------------------                          *
'*     1.000000         8.999984                           *
'*     2.000000        80.999881                           *
'*     3.000000       360.999496                           *
'*     4.000000      1088.998512                           *
'*     5.000000      2600.996484                           *
'*     6.000000      5328.992838                           *
'*     7.000000      9800.986875                           *
'*     8.000000     16640.977767                           *
'*     9.000000     26568.964559                           *
'*    10.000000     40400.946171                           *
'*  -----------------------------                          *
'*                                                         *
'}**********************************************************
'{ PROGRAM EQUDIFPC
WindowStyle 24:Font 2:randomize:Window 0,0-%maxx,%maxy-40
declare tx![100],ty![100],x!,xi!,xf!,yi!,  fi&,n&,i&
if 1   '  0=Manual input, 1=Example
  xi!=0  :print "\n\n X_begin: ";xi!
  xf!=10 :print " X_end  : ";xf!
  yi!=1  :Print " Y_begin: ";yi!
  n& =10 :print " Number of points inbetween: ";n&
  fi&=10 :print " Granularity between points: ";fi&
else
  print "\n\n Input x begin: ";  :input xi!
  print " Input x end  : ";      :input xf!
  Print " Input y begin: ";      :input yi!
  print " Input number of points: ";:input n&
  print " Input Granularity: ";  :input fi&
endif


equadiff_pc tx![],ty![],xi!,xf!,yi!,n&,fi&


print ' Results:
print "      X                Y      "
print " -----------------------------"
WhileLoop n&:i&=&Loop 
  print format$("%g",tx![i&]),format$("%g",ty![i&])
endwhile
print " ----------------------------- "
waitinput
END 
'} End of main program


' Hier die Definition der Differentialgleichung f: 
' f: y'=( 4xy + 4x*sqrt(y))/(1+x^2)
proc f :parameters x!,y! 
declare f!
if y!>=0


f! = ( 4*x!*y!+4*x!*sqrt(y!) ) / ( 1+sqr(x!) )
endif
return f!
endproc


proc runge_kutta :parameters h!,x!,y!
' classical Runge-Kutta method of order 4
declare a!,b!,c!,d!,f!
  a!=h!*f(x!,y!)
  b!=h!*f(x!+h!/2, y!+a!/2)
  c!=h!*f(x!+h!/2, y!+b!/2)
  d!=h!*f(x!+h!,y!+c!)
y!=y!+(a!+2*b!+2*c!+d!)/6
z!=y!
endproc


'**********************************************************
'{             Prediction-correction method               *
'* ------------------------------------------------------ *
'* INPUTS:   xi      begin x value                        *
'*	         xf      end x value                          *
'*           y1      begin y value (at x=xi)              *
'*           m       number of points to calculate        *
'*           fi      finesse (number of intermediate      *
'*                   points (for example 10)              *
'* OUTPUTS:                                               *
'*           tx      table of x values (m values)         *
'*           ty      table of y values (m values)         *
'*                                                        *
'* DESCRIPTION:                                           *
'* The algorithm has the following steps:                 *
'*     1. calculate y2, y3, y4 using a classical Runge-   *
'*        Kutta method of order 4.                        *
'*     2. from point (x4, y4), first estimate y(n+1) by   *
'*        PREDICTION formula:                             *
'*  y(n+1)=y(n)+h/24(55y'(n)-59y'(n-1)+37y'(n-2)-9y'(n-3) *
'*                                                        *
'*  then continue with CORREKTION formula:                *
'*  y(n+1)=y(n) + h/24(9y'(n+1)+19y'(n)-5y'(n-1)+y'(n-2)  *
'*                                                        *
'*  substituting the y'(n+1) = f(x(n+1),y(n+1)) by the    *
'*        estimated value of y(n+1), until convergence    *
'*        is obtained.                                    *
'}*********************************************************


Proc equadiff_pc :parameters tx![],ty![],xi!,xf!,yi!,m&,fi&
declare ni&,h!,w!,x!,y!,z!,p![3] , k&,i&,j&


z!=yi! : case (m&>100) or (fi&<1) : Goto "retour"


h!=(xf!-xi!)/fi&/m&
p![3]=f(xi!,yi!) : tx![0]=xi! : ty![0]=yi!


k&=0
whileloop m&:i&=&Loop : ni&=(i&-1)*fi&-1
  whileloop fi&:j&=&Loop : x!=xi!+h!*(ni&+j&) : inc k&
    if k&<4
      runge_kutta h!,x!,z!
      x!=x!+h!  '<<< wird in Fortran90 in runge_kutta selbst als Rückparameter erledigt! 
      p![3-k&]=f(x!,z!)
    else
      x!=x!+h!
      w!=z! +h!/24*(55*p![0]-59*p![1] +37*p![2]  -9*p![3])
' y(n+1)=y(n)+h/24*(55*y'(n)-59y'(n-1)+37*y'(n-2)-9*y'(n-3)
      Repeat
        y!=w!
        w!=z!+h!/24*(9*f(x!,y!)+19*p![0]-5*p![1] + p![2])
'   y(n+1)=y(n)+h/24(9*y'(n+1) +19*y'(n)-5*y'(n-1)+y'(n-2)
      Until abs(y!-w!)<val("1e-10")
      z!=w! : p![3]=p![2] : p![2]=p![1] : p![1]=p![0]
      p![0]=f(x!,z!)
    endif
  endwhile 'j loop
  tx![i&]=x!
  ty![i&]=z!
endwhile 'i loop
endproc

Proc DEQHelp
font 2:Cls:Print "\n\n============================================================"
Print "                        D E Q H e l p"
Print "============================================================"
Print "\n CLEAR q&[] ' natives XProfan ":Print " sz&  = SizeOf(q&[]) ' sz=5 bedeutet: 0,1,2,3,4 "
Print " SetSize q&[],sizeof(q&[])+1 ' erweitert Array um 1 Stelle "
Print " MakeIndexDEQ(q&[],n&) ":Print " MakeRevertIndexDEQ(q&[],n&) "
Print " MakeConstDEQ(q&[],n&,wert&) ":Print "\n ShowDEQ q&[],nw&   ' Temporary NumWidth "
Print "\n Wert&= PopBtm(q&[]) ' -999999999 = OutOfRange Error "
Print " Wert&= PopTop(q&[]) ' or -999999999 ":Print " PutTop Wert&,q&[] "
Print " PutBtm Wert&,q&[] ":Print "\n Wert&= CopyBtm(q&[]) "
Print " Wert&= CopyTop(q&[]) ' -999999999 ":Print " Wert&= CopyFromDEQ(q&[],n&) "
Print "\n LogVgl&=CmpInDEQ(q&[],n&,m&)' n=Pos1,m=Pos2: <1,'=0',>-1 "
Print "\n NisM&= SwapInDEQ(q&[],n&,m&)' (N=M)=1, -999999999 "
Print " Wert&= DelFromDEQ(q&[],n&) ":Print " DEQHelp ' diesen Text ausgeben":waitinput :endproc
:proc MakeIndexDEQ :parameters q&[],n&:SetSize q&[],n&:q&[]=&index+1:endproc
:proc MakeRevertIndexDEQ :parameters q&[],n&:SetSize q&[],n&:q&[]=n&-&index:endproc
:proc MakeConstDEQ :parameters q&[],n&,wert&:SetSize q&[],n&:q&[]=Wert&:endproc
:Proc PopBtm :parameters q&[]:var sz&=sizeof(q&[]):if sz&>0:var wert&=q&[0]
q&[]=q&[&index+(&index<(sz&-1))]:setsize q&[],sizeof(q&[])-(sz&>0):return wert&
else :return int(-999999999):endif :endproc
:proc PopTop :parameters q&[]:declare wert&:var sz&=sizeof(q&[]):if sz&>0:wert&=q&[sz&-1]
else :wert&=int(-999999999):endif :setsize q&[],sz&-(sz&>0):return wert&:endproc
:Proc PutTop :parameters Wert&,q&[]:var sz&=sizeof(q&[]):inc sz&:setsize q&[],sz&:q&[sz&-1]=Wert&:endproc
:Proc PutBtm :parameters Wert&,q&[]:declare tmp&[]:var sz&=sizeof(q&[]):inc sz&
setsize q&[],sz&:tmp&[]=q&[]:q&[]=tmp&[&index-(&index>0)]:clear tmp&[]:q&[0]=Wert&:endproc
:Proc CopyBtm :parameters q&[]:var sz&=sizeof(q&[]):if sz&>0:return q&[0]
else :return int(-999999999):endif :endproc
:proc CopyTop :parameters q&[]:declare wert&:var sz&=sizeof(q&[])
if sz&>0:wert&=q&[sz&-1]:else :wert&=int(-999999999):endif :return wert&:endproc
:proc ShowDEQ :parameters q&[],nw&:var tmp&=get("numwidth"):set("numwidth",nw&)
var sz&=sizeof(q&[]):print tab(2);"n=";sz&;":  ";:whileloop 0,sz&-1:print q&[&Loop],
endwhile :print:set("numwidth",tmp&):endproc
:proc CopyFromDEQ :parameters q&[],n&:declare wert&:var sz&=sizeof(q&[])
if (sz&>0) and (n&>=0) and (n&<sz&):wert&=q&[n&]:else :wert&=int(-999999999):endif :return wert&:endproc
:Proc CmpInDEQ :parameters q&[],n&,m&
var sz&=sizeof(q&[]):if (sz&>0) and (n&>=0) and (n&<sz&) and (m&>=0) and (m&<sz&)
return int((q&[n&]<q&[m&])-(q&[n&]>q&[m&])):else :return int(-999999999):endif :Endproc
:Proc SwapInDEQ :parameters q&[],n&,m& 'Swaps two elements in Queue
var sz&=sizeof(q&[]):if (sz&>0) and (n&>=0) and (n&<sz&) and (m&>=0) and (m&<sz&):var tmp&=q&[n&]
q&[n&]=q&[m&]:q&[m&]=tmp&:return int(n&<>m&):else :return int(-999999999):endif :endproc
:proc DelFromDEQ :parameters q&[],n&:var sz&=sizeof(q&[]) :declare wert&
if (sz&>0) and (n&>0) and (n&<sz&):wert&=q&[n&]:q&[]=q&[&index+(&index>=n&)*(&index<(sz&-1))]
dec sz&:setsize q&[],sz&:else :wert&=int(-999999999):endif :return wert&
endproc

WindowTitle "Fixed-StringArray Memory Testfacility"
Cls:font 0:Randomize:set("randseed",1)
Declare a%,a&,b%[],a$[9],flg%
a$[]=left$(\
chr$(56+rnd(32))+chr$(56+rnd(32))+chr$(56+rnd(32))+chr$(56+rnd(32))\
+chr$(56+rnd(32))+chr$(56+rnd(32))+chr$(56+rnd(32))+chr$(56+rnd(32))\
+chr$(56+rnd(32))+chr$(56+rnd(32))+chr$(56+rnd(32))+chr$(56+rnd(32))\
+chr$(56+rnd(32))+chr$(56+rnd(32)),3+&index)
proc prt
case flg%:print tab(5);&index;".",a$[&index]
inc flg%:return a$[&index]
endproc
flg%=0:a$[]=prt()
flg%=0:a$[]=prt()
a%=a$[]:font 2:print " Das an",a%;". Stelle deklarierte Array!":font 0


whileloop 0,9
a%=addr(a$[&Loop]):print "   ";a%,"len:";long(a%,-4),string$(a%,0),"Typ:";long(a%,-8),
print "Endbytes:",byte(addr(a$[&Loop])+long(a%,-4),0);",";byte(addr(a$[&Loop])+long(a%,-4)+1,0);
print ",";byte(addr(a$[&Loop])+long(a%,-4)+2,0);",";byte(addr(a$[&Loop])+long(a%,-4)+3,0)
endwhile


waitinput
end

WindowTitle "VELOCITY VERLET DEMO am Beispiel Lineares Feder-Masse-Dämpfungssystem"
'(CL) CopyLeft 2014-08 by P.Specht, Wien
'Q: http://en.wikipedia.org/wiki/Verlet_integration#Velocity_Verlet
WindowStyle 24:Window 0,0-%maxx,%maxy
var xh&=width(%hwnd)\2:var yh&=height(%hwnd)\2
Declare m!,dt!,x1!,x2!,v1!,v15!,v2!,a1!,a2!


proc A :parameters x!,v!
var l!=200 ' Federlänge m/1000
var fk!=25 ' Federkonstante N/m 'Kraft pro Auslenkung
var rk!=0.2' Reibungskonstante N*s/m
return (l!-x!)*fk!/m! - v!*rk!
'Beschleunigung = Kraft / Masse - Reibungsdämpfung
endproc


proc DRAW
usepen 0,(xh&-x2!)/30  ,255
cls
line xh&-500,yh& - xh&+x2!-200,yh&
usepen 0,37,rgb(0,0,255)
line xh&+x2!-200,yh& - xh&+x2!-199,yh&
endproc


'Startwerte:
dt!=0.03 's
m!=1     'kg


'Einfachstes Zeichenmodell (mit Cls, ohne MemCopy)


'Hauptschleife
Repeat :waitinput 33




'VELOCITY VERLET (aka Leapfrog)
v15!= v1!+a1!/2*dt!
x2! = x1!+v15!*dt!
a2! = A(x2!,v15!)
v2! = v15!+dt!*a2!/2
DRAW
a1!=a2!:v1!=v2!:x1!=x2! ' Umgreifen für nächsten Schritt




until %key=27
END

WindowTitle mkstr$(" ",26)+upper$("AKIMA-Interpolation mittels Kubischer Splines [H. Akima, 1978]")
'  Portiert nach XProfan11 (CL) CopyLeft 2014-08 by P. Specht, Vienna (Austria).
'{ Keine wie auch immer geartete Garantie. Nutzung auf alleinige Gefahr des Anwenders!
'
' Übersetzungsbasis: Turbo Pascal Handbuch 1974; Für Detailklärungen wurde verwendet:
' F90/95: https://www.tacc.utexas.edu/documents/13601/162125/fortran_class.pdf
' F95: http://www.mrao.cam.ac.uk/~rachael/compphys/SelfStudyF95.pdf
'
' Pascal-DemoSource: http://jean-pierre.moreau.pagesperso-orange.fr/Pascal/akima_pas.txt
' Orig. 2D-Algorithmus: Hiroshi Akima,U.S. Department of Commerce, NTIA/ITS, F90-Version of 1995/08, 
' veröffentlicht in "Hiroshi Akima, 'A Method of Bivariate Interpolation and Smooth Surface Fitting
' for Irregularly Distributed Data Points", ACM Transactions on Math. Software Vol.4 No.2, June 1978
' pp. 148-159. Copyright 1978, Association for Computing Machinery, Inc.
' Fortran-90: TOMS 22.yr,No3 (Sept. 1996) p.357ff,  http://calgo.acm.org/760.gz (Abfrg. 2014-08)
' Test Data: http://cran.r-project.org/web/packages/akima/akima.pdf
'
'********************************************************
'*   Program to demonstrate the Akima spline fitting    *
'*       of Function SIN(X) in double precision         *
'* ---------------------------------------------------- *
'*   Reference: BASIC Scientific Subroutines, Vol. II   *
'*   By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].  *
'*                                                      *
'*             Pascal Version by J.-P. Moreau, Paris.   *
'*                       (www.jpmoreau.fr)              *
'* ---------------------------------------------------- *
'* SAMPLE RUN:                                          *
'*                                                      *
'* Akima spline fitting of SIN(X):                      *
'*                                                      *
'*   X   SIN(X) HANDBOOK   AKIMA INTERPOLATION   ERROR  *
'* ---------------------------------------------------- *
'* 0.00    0.0000000          0.0000000       0.0000000 *
'* 0.05    0.0499792          0.0500402      -0.0000610 *
'* 0.10    0.0998334          0.0998435      -0.0000101 *
'* 0.15    0.1494381          0.1494310       0.0000072 *
'* 0.20    0.1986693          0.1986459       0.0000235 *
'* 0.25    0.2474040          0.2474157      -0.0000118 *
'* 0.30    0.2955202          0.2955218      -0.0000016 *
'* 0.35    0.3428978          0.3428916       0.0000062 *
'* 0.40    0.3894183          0.3894265      -0.0000081 *
'* 0.45    0.4349655          0.4349655       0.0000000 *
'* 0.50    0.4794255          0.4794204       0.0000051 *
'* 0.55    0.5226872          0.5226894      -0.0000021 *
'* 0.60    0.5646425          0.5646493      -0.0000068 *
'* 0.65    0.6051864          0.6051821       0.0000043 *
'* 0.70    0.6442177          0.6442141       0.0000035 *
'* 0.75    0.6816388          0.6816405      -0.0000017 *
'* ---------------------------------------------------- *
'*                                                      *
'*******************************************************}
'}


PROC Interpol_Akima
'********************************************************
'*       Akima Spline Fitting Subroutine                *
'* -----------------------------------------------------*  
'* The input table is X(i), Y(i), where Y(i) is the     *
'* dependant variable. The interpolation point is xx,   *
'* which is assumed to be in the interval of the table  *
'* with at least one table value to the left, and three *
'* to the right. The interpolated returned value is yy. *
'* n is returned as an error check (n=0 implies error). *
'* It is also assumed that the X(i) are in ascending    *
'* order!                                               *
'********************************************************
parameters i&
n&=1
' {special case xx=0}
if xx!=0: yy!=0:return:endif
' {Check to see if interpolation point is correct}
if (xx!<X![1]) or (xx!>=X![iv&-3]): n&=0:return:endif
X![0]=2*X![1]-X![2]
' {Calculate Akima coefficients, a and b}
whileloop iv&-1:i&=&Loop
'   {Shift i to i+2}
  XM![i&+2]=(Y![i&+1]-Y![i&])/(X![i&+1]-X![i&])
endwhile
XM![iv&+2]=2*XM![iv&+1]-XM![iv&]
XM![iv&+3]=2*XM![iv&+2]-XM![iv&+1]
XM![2]=2*XM![3]-XM![4]
XM![1]=2*XM![2]-XM![3]
whileloop iv&:i&=&Loop
  a!=ABS(XM![i&+3]-XM![i&+2])
  b!=ABS(XM![i&+1]-XM![i&])
  if (a!+b!)<>0
    Z![i&]=(a!*XM![i&+1]+b!*XM![i&+2])/(a!+b!)
  else
    Z![i&]=(XM![i&+2]+XM![i&+1])/2
  endif
endwhile   
' {Find relevant table interval}
i&=0
repeat
  inc i&
until xx!<=X![i&]
dec i&
' {Begin interpolation}
b!=X![i&+1]-X![i&]
a!=xx!-X![i&]
yy!=Y![i&]+Z![i&]*a!+(3*XM![i&+2]-2*Z![i&]-Z![i&+1])*a!*a!/b!
yy!=yy!+(Z![i&]+Z![i&+1]-2*XM![i&+2])*a!*a!*a!/(b!*b!)
EndProc


'{ Main Program AKIMA-SINTEST
Windowstyle 24:font 0:set("decimals",8)
CLS 'Window 0,0-%maxx,%maxy
var fmt$=" 0.0000000;-0.0000000"


var SIZE& = 14


declare X![SIZE&+1],Y![SIZE&+1],XM![SIZE&+4],Z![SIZE&+1]
declare a!,b!,xx!,yy!, i&,iv&,n&


iv&=14  '{Number of points in table}
CLS
Font 2:Print "\n";tab(10);"  Hiroshi Akima's Interpolation mit kubischen Splines"
Font 0:print tab(10);" "+mkstr$("-",53)
print tab(10);" Basis sind 14 unregelmäßig verteilte Stützwerte einer "
print tab(10);" Eigenbau-Sinusfunktion. Errechnet wird der absolute   "
Print tab(10);" Fehler (als Differenz zur internen Sin(x)-Funktion).  "
print tab(10);" Im Programmtext sind die Vergleichswerte zu finden!   \n"
' -----------------------------------------------------------------
' Input y=Sine(x) table
' -----------------------------------------------------------------
' Sine table values from  Handbook of mathematical functions
' by M. Abramowitz and I.A. Stegun, NBS, june 1964
' -----------------------------------------------------------------
  X![1]=0.000:Y![1]=0.00000000  :  X![2]=0.125:Y![2]=0.12467473
  X![3]=0.217:Y![3]=0.21530095  :  X![4]=0.299:Y![4]=0.29456472
  X![5]=0.376:Y![5]=0.36720285  :  X![6]=0.450:Y![6]=0.43496553
  X![7]=0.520:Y![7]=0.49688014  :  X![8]=0.589:Y![8]=0.55552980
  X![9]=0.656:Y![9]=0.60995199  :  X![10]=0.721:Y![10]=0.66013615
  X![11]=0.7853981634           : Y![11]=0.7071067812
  X![12]=0.849:Y![12]=0.75062005:  X![13]=0.911:Y![13]=0.79011709
  X![14]=0.972:Y![14]=0.82601466
' -----------------------------------------------------------------


print "\n";mkstr$("-",77)
font 2
print tab(3+3);"X";tab(%pos+15);"SIN(X)";\
        tab(%pos+3);"AKIMA-Stützwertinterpolation";tab(%pos+2);"ABS.FEHLER"
font 0
print mkstr$("=",77)


' {main loop}
xx!=0
whileloop 16:i&=&Loop
  Interpol_Akima(i&)
  print tab(4);xx!;tab(%pos+6);SIN(xx!);tab(%pos+9);yy!;tab(%pos+11);format$(fmt$,SIN(xx!) - yy!)
  casenot i& mod 4:print
  xx! = xx! + 0.05
endwhile
' {print footer}
print mkstr$("-",77)
beep
waitinput
END 
'} of file akima.prf