#Crib: http://www.math.uni-wuppertal.de/~xsc/literatur/PXSCENGL.pdf # INT n = 4; # x, xd, y, yd # []STRING cmm = ("x ", "xd ", "y ", "yd "); MODE R = LONG LONG REAL; MODE RF = PROC(R)R; RF rsqrt = long long sqrt; MODE SV = [n]R; PROC hlf = (SV x)SV: (SV a; FOR j TO n DO a[j] := 0.5*x[j] OD; a); PROC prs = (STRING m, SV x)VOID: (print((m, newline)); FOR j TO n DO print((cmm[j], x[j], newline)) OD); OP * = (R s, SV v)SV: (SV a; FOR j TO n DO a[j] := s*v[j] OD; a); OP * = (REAL s, SV v)SV: (SV a; FOR j TO n DO a[j] := s*v[j] OD; a); OP + = (SV x, y)SV: (SV a; FOR j TO n DO a[j] := x[j] + y[j] OD; a); # This code is to solve the ordinary differential equation dy/dx = f(x, y) # # in particular a Newtonian prabolic orbit # PROC f = (R tau, SV y)SV: (R d2 = y[1]*y[1] + y[3]*y[3], d = rsqrt(d2), dm3 = 1.0/(d*d2); (SV a; FOR i BY 2 TO n DO a[i] := y[i+1]; a[i+1] := -y[i]*dm3 OD; a)); PROC step = (R dx, x, SV y)SV: ( SV k1 = dx*f(x,y); SV k2 = dx*f(x+dx/2, y + hlf(k1)); SV k3 = dx*f(x+dx/2, y + hlf(k2)); SV k4 = dx*f(x+dx, y+k3); y + 1.0/R(6)*(k1+2.0*k2 + 2.0*k3 + k4)); (SV x := (1, 0, 0, rsqrt(2)); PROC rep = (INT j)VOID: (print((j, newline)); prs("dvs", x); prs("psh", f(0, x)); # This orbit should be 4x + y^2 = 4 # print(("parabola ", 4.0*x[1] + x[3]*x[3] - 4.0, newline))); INT c = 10000, top= 55000; R dx = 1/c; rep(0); FOR k TO top DO x := step (dx, 0, x); (k=1 OR (k MOD c) = 0 OR k=top | rep(k)) OD) # +0 dvs x +1.00000000000000000000000000000000000000000000000000000000000000e +0 xd +0.00000000000000000000000000000000000000000000000000000000000000e +0 y +0.00000000000000000000000000000000000000000000000000000000000000e +0 yd +1.41421356237309504880168872420969807856967187537694807317667974e +0 psh x +0.00000000000000000000000000000000000000000000000000000000000000e +0 xd -1.00000000000000000000000000000000000000000000000000000000000000e +0 y +1.41421356237309504880168872420969807856967187537694807317667974e +0 yd +0.00000000000000000000000000000000000000000000000000000000000000e +0 parabola +0.00000000000000000000000000000000000000000000000000000000000000e +0 +1 dvs x +9.99999995000000016666666604166667083333330647786472737630132853e -1 xd -9.99999993333333395833332588541675367838442325847475158678922501e -5 y +1.41421356001607245368536503017788742465402489787791994205820225e -4 yd +1.41421355530202729586177791616761541682110383262411492058686998e +0 psh x -9.99999993333333395833332588541675367838442325847475158678922501e -5 xd -9.99999980000000291666662833333381249999418652350667046359779982e -1 y +1.41421355530202729586177791616761541682110383262411492058686998e +0 yd -1.41421353880286933628698666311526347883250219011089942454470199e -4 parabola +5.55555555555555567925347031250001460780000000000000000000000000e -26 +10000 dvs x +6.08721781282468751692660505803174200265888275298245741252329202e -1 xd -6.35834147689268604754633457052697722842831821201157217756724281e -1 y +1.25104471337763343240798572367813295315942208423974820689400843e +0 yd +1.01648508784727860505754688668825036381473443944872199550603365e +0 psh x -6.35834147689268604754633457052697722842831821201157217756724281e -1 xd -2.26035620333335739558702594732548478187805279378849917518198682e -1 y +1.01648508784727860505754688668825036381473443944872199550603365e +0 yd -4.64548298661639649391653880805302416042850401009555680484639839e -1 parabola -5.95464417034105218716849723792691633217301774400000000000000000e -18 +20000 dvs x -8.08594603928763757933625488764077898515421823920041196779860918e -2 xd -7.06572714825347830630488227851612025471746262620298675203375059e -1 y +2.07928782076255750847655386436364664839057631788351550655651463e +0 yd +6.79629542163354297188820805010140567709913083027041187376821185e -1 psh x -7.06572714825347830630488227851612025471746262620298675203375059e -1 xd +8.97434148374444182264227070510602640848572109056195159360225951e -3 y +6.79629542163354297188820805010140567709913083027041187376821185e -1 yd -2.30773725867677725146277617755425765513526988514796059616794295e -1 parabola -2.37422658302811365871405444748382674333110275880000000000000000e -17 +30000 dvs x -7.75726623466793165272655302622630841738590072468473895073458386e -1 xd -6.78932126976413522689937460140925159651921588943738147538280182e -1 y +2.66512785694554860949873044103396424209322859836076105579762108e +0 yd +5.09493100083029014757365854654172017849579191784035050717589420e -1 psh x -6.78932126976413522689937460140925159651921588943738147538280182e -1 xd +3.62725947696178928252981859682932813239334441397156109197553253e -2 y +5.09493100083029014757365854654172017849579191784035050717589420e -1 yd -1.24620065677485785192611698602253548938331536247330496909699995e -1 parabola -4.74422936474959736487132028423917795430874115820000000000000000e -17 +40000 dvs x -1.43628491597273211684678882829342515976918628469014727463476569e +0 xd -6.42376833488174989335804721722704890545573918338482975748293954e -1 y +3.12172062553504752155487112627847502426474078163447188499339863e +0 yd +4.11553057140130699773201921208241375985543944699531326779816550e -1 psh x -6.42376833488174989335804721722704890545573918338482975748293954e -1 xd +3.53975417216011261571747975754473095748165772150660257882349804e -2 y +4.11553057140130699773201921208241375985543944699531326779816550e -1 yd -7.69354567862477486057116796756975410235194880031846863823745653e -2 parabola -7.53149329346107004030000776205826127273329123430000000000000000e -17 +50000 dvs x -2.06170354394960126602381594313020227691951010114206100324377622e +0 xd -6.09239908725110659833712053086914023276163104929951802120670241e -1 y +3.49954485266275805132799629698339204891903173440772215224831487e +0 yd +3.48182369065250270162744478431629570663697716504174329694254389e -1 psh x -6.09239908725110659833712053086914023276163104929951802120670241e -1 xd +3.07681618748136252044925044702652318158812779206652061878570544e -2 y +3.48182369065250270162744478431629570663697716504174329694254389e -1 yd -5.22260161170536662975510253291276874320869652500730161475779558e -2 parabola -1.06363133128754140584937939304370127584785134320000000000000000e -16 +55000 dvs x -2.36257531466267626833087516671515017160637782612458094710346185e +0 xd -5.94440061483565566965427277705964364264949417541064577397498494e -1 y +3.66746523618843658568955312785486728890940354584215292865142018e +0 yd +3.24169432128748264472828432904265944027938932771592227781244397e -1 psh x -5.94440061483565566965427277705964364264949417541064577397498494e -1 xd +2.84548859403596451176808658424462793898161407594436910605428590e -2 y +3.24169432128748264472828432904265944027938932771592227781244397e -1 yd -4.41709960898647730905265350913433630999084398489011679767589481e -2 parabola -1.22874049950468686316834336937520843658627949764000000000000000e -16 #