rem Circular Surface Waves .... Rev 9.0 rem (Assembly Language version) rem A J Tooth / 20th August 2003 rem Modified with Assembly Routines April 2005 *FLOAT 64 himem=lomem + 10000000 rem Set up use of Full Screen proc_fullscreen(22) rem Parameter Setup proc_param print " Calculating Bessel Function Array " proc_ArrBess(Xmax,K) cls rem Read in Weights proc_wghtcalc rem Further Setup origin (Xlim%/2),0 |dk%=(K-1)/qt% rem Time Loops for ?tim%=0 to (?stp%-1) Tm=(1+?tim%)*Scal%/?stp% rem Calculate Waveforms print tab(5,5);" Calculating Wave for time-slot ";1+?tim%;"/";?stp% proc_calc next print tab(5,10);"Calculating. Please wait a few seconds...." *REFRESH *REFRESH OFF rem Display Resultant Surface Wave proc_disp *REFRESH ON a$=get$ : quit end rem End of Program ....................................................... rem----------------------------------------------------------------------- rem Traps Computational Range Errors ONLY 450 if err=21 or err=23 or err=18 or err=20 or err=22 or err=24 then goto 1560 else error endif goto 1620 rem----------------------------------------------------------------------- rem Parameter Setup def proc_param Xlim%=2048 : Ylim%=1536 rem Enter Parameters input " How many calculated weights (24)",qt% if qt%=0 then qt%=24 input " Enter calculated weight factor (.01)",c if c=0 then c=0.01 input " Timescale (200) ",Scal% if Scal%=0 then Scal%=200 input " Number of time-steps (50) ",sp% if sp%=0 then sp%=50 if sp%>255 then sp%=255 input " Enter largest x-value (75) ",Xmax if Xmax=0 then Xmax=75 input " Enter max spread of wave-number,K (1.5) [Must be >1] ",K if K<1.0 then K=1.5 K*=1.0 input " Wave Height Factor (500)",Wh% if Wh%<=0 then Wh%=500 input " Enter height (300)",vh% if vh%<=0 then vh%=300 input " Enter distance away (400)",hd% if hd%<=0 then hd%=400 rem Dimension various Arrays dim dm% 3, dml% 3, arrsize% 3, ref% 3, stp% 0 dim vht% 3, hdis% 3, Wht% 3, dk% 7, ptts% 7 !vht%=vh% : !hdis%=hd% ?stp%=sp% !dm%=Xlim%/16 !dml%=!dm%-1 !arrsize%=16*(!dm%*!dm%) !Wht%=Wh% dim pnts(!dm%,?stp%), scl% 3, x% 3, y% 3, vy% 3, col% 3, cd% 0 dim pmts% (8*!dm%*?stp%), sef% 3, tef% 3 dim alpbet% !arrsize% dim Fniw% (!arrsize%/2) dim prs% !arrsize%, Sf% 7, strt% 3, fset% 3 dim bess(10000) dim Wght%(qt%), Tot% 3, j% 0 !scl%=4000 dim drw% 500, rr% 3, cc% 3, px% 7, py% 7, fw% 7 dim Fdrw% 300, AX% 7, AY% 7 dim cent% 3, q% 3, q2% 7, t% 0, ang% 3, ex% 3, ey% 3 dim circ% 300, norm% 100, calc% 300, tim% 0 dim Mdrw% 200, alpha% 400, rat% 7, beta% 400 dim xd% 7, yd% 7, zd% 7, Ans% 7, dmp% 7, Rr% 7 rem Set up Arcos array proc_arcos_setup rem Dual-pass assembly, in case of labels for pass%=0 to 2 step 2 proc_Mdrw(pass%) proc_drw(pass%) proc_circ(pass%) proc_Fdrw(pass%) proc_norm(pass%) proc_calc(pass%) proc_alpha(pass%) proc_beta(pass%) next pass% cls endproc rem----------------------------------------------------------------------- rem Display Waves def proc_disp !cent%=!dm%/2 for ?t%=0 to (?stp%-1) rem Distribute the wave shape in a circular fashion call circ% rem Calculate perspective proc_pers rem Display the wave in perspective cls sys "GetCurrentProcess" to hprocess% sys "SetPriorityClass", hprocess%, &100 rem Displays Terrain call Mdrw% sys "GetCurrentProcess" to hprocess% sys "SetPriorityClass", hprocess%, &20 *REFRESH next endproc rem----------------------------------------------------------------------- rem Calculate Wave Points def proc_calc for !x%=0 to (!dm%-1) xac=!x%*Xmax/!dm% |ptts%=0.0 rem Combine Fourier Integral Elements for ?j%=0 to qt% on error local goto 450 1560 k=(1+(?j%*|dk%))/2 w=sqr(10*k) xpass=k*xac |ptts%+=Wght%(?j%)*w*(sin(w*Tm))*fn_GetBess(xpass) next call calc% 1620 next endproc rem----------------------------------------------------------------------- rem Calculated Weights Procedure. Uses exp(-x²) function. def proc_wghtcalc local p%,v !Tot%=0 for p%=0 to qt% v=p%-(qt%/2) Wght%(p%)=int(1000*(exp(-c*(v*v)))) !Tot%+=Wght%(p%) next p% endproc rem----------------------------------------------------------------------- rem Set up use of Full Screen def proc_fullscreen(N%) sys "GetSystemMetrics", 0 to xscreen% sys "GetSystemMetrics", 1 to yscreen% sys "SetWindowLong",@hwnd%,-16,&16000000 sys "SetWindowPos",@hwnd%,-1,0,0,xscreen%,yscreen%,0 mode N% mouse off : off : rem Turns off the Mouse Pointer and the Cursor endproc rem----------------------------------------------------------------------- rem Traps Computational Range Errors ONLY def proc_error if err=21 or err=23 or err=18 or err=20 or err=22 or err=24 then rem Ignore in this case skp%=1 else skp%=0 : restore error endif endproc rem----------------------------------------------------------------------- rem Calculate Pre-set Reference Array for Bessel Function def proc_ArrBess(Xmax,K) local a%,x for a%=0 to 10000 x=a%*(K+1)*Xmax/10000 bess(a%)=fn_Bess(x) next a% endproc rem----------------------------------------------------------------------- rem Calculates Bessel Function Values def fn_Bess(X) local dTh,Cumm,Th,t% Sc%=100 dTh=pi/Sc% Cumm=0 for t%=1 to Sc% Th=pi*(t%-0.5)/Sc% Cumm=Cumm + dTh*cos(X*sin(Th)) next t% Cumm=Cumm/pi =Cumm rem----------------------------------------------------------------------- rem Gets the appropriate Bessel Function value from the pre-set array def fn_GetBess(X) local Ref% Ref%=int(10000*X/((K+1)*Xmax)) =bess(Ref%) rem----------------------------------------------------------------------- rem Calculates perspective def proc_pers sys "GetCurrentProcess" to hprocess% sys "SetPriorityClass", hprocess%, &100 |xd%=1.0*!vht% : |yd%=1.0*!hdis% : |zd%=0.0 call norm% |Rr%=|Ans% for !y%=0 to (!dm%-1) for !x%=0 to (!dm%-1) skp%=0 on error local proc_error if skp%=0 then call alpha% call beta% endif next next sys "GetCurrentProcess" to hprocess% sys "SetPriorityClass", hprocess%, &20 endproc rem----------------------------------------------------------------------- rem Assembly Language Routine for Plotting def proc_drw(opt%) P%=drw% [opt opt% finit mov eax,[ref%] mov ebx,prs% add ebx,eax add ebx,8 fld qword [ebx] sub ebx,8 fld qword [ebx] fild dword [scl%] fmulp st1,st0 fistp dword [rr%] fild dword [scl%] fmulp st1,st0 fistp dword [cc%] mov al,25 \Plot XX,rr%,cc% call "oswrch" mov al,[cd%] call "oswrch" mov eax,[rr%] shl eax,1 mov bx,ax mov al,bl call "oswrch" mov al,bh call "oswrch" mov eax,[cc%] shl eax,1 mov bx,ax mov al,bl call "oswrch" mov al,bh call "oswrch" \Plot XX,rr%,cc% ret ] endproc rem----------------------------------------------------------------- rem Assembly Language Routine for Drawing def proc_Fdrw(opt%) P%=Fdrw% [opt opt% mov edx,0 mov [x%],edx \Initialise x-loop counter .loopx mov eax,[vy%] imul eax,[dm%] add eax,[x%] shl eax,3 mov ebx,Fniw% add ebx,eax finit fld qword [Sf%] fld st0 fld1 fld1 faddp st1,st0 fdivp st1,st0 fld qword [ebx] faddp st1,st0 fdivrp st1,st0 fild dword [tff%] fmulp st1,st0 fistp dword [col%] mov al,19 \Changes the rgb vlue f0r colr 10 call "oswrch" mov al,10 call "oswrch" mov al,16 call "oswrch" mov al,0 call "oswrch" mov al,0 call "oswrch" mov al,[col%] call "oswrch" mov al,18 \Calls gc0l routine call "oswrch" mov al,0 call "oswrch" mov al,10 call "oswrch" mov eax,[vy%] imul eax,[dm%] add eax,[x%] shl eax,4 mov [ref%],eax call drw% dec dword [vy%] mov eax,[vy%] imul eax,[dm%] add eax,[x%] shl eax,4 mov [ref%],eax call drw% inc dword [vy%] inc dword [x%] mov edx,[x%] cmp edx,[dml%] jle near loopx jmp over .tff% DD 255 .over ret ] endproc rem ============================================================ rem Assembly Language Routine for distributing waveform rem cylindrically def proc_circ(opt%) P%=circ% [opt opt% mov edx,0 mov [q%],edx \Initialise q-loop counter .loopq mov eax,0 \Calculation 0f sef% reference t0 pnts array mov al,[t%] imul eax,[dm%] add eax,[q%] shl eax,3 mov ebx,pmts% add ebx,eax mov [sef%],ebx finit fild dword [q%] fld1 fld1 faddp st1,st0 fdivp st1,st0 fstp qword [q2%] mov edx,0 mov [ang%],edx \Initialise ang-loop counter .loopang fldpi fild dword [ang%] fmulp st1,st0 fild dword [fivhun%] fdivp st1,st0 fld st0 fcos fld qword [q2%] fmulp st1,st0 fild dword [cent%] faddp st1,st0 fistp dword [ex%] fsin fld qword [q2%] fmulp st1,st0 fild dword [cent%] faddp st1,st0 fistp dword [ey%] mov eax,[ey%] imul eax,[dm%] add eax,[ex%] shl eax,3 mov [tef%],eax mov esi,[sef%] fld qword [esi] mov ebx,Fniw% add ebx,eax fstp qword [ebx] inc dword [ang%] mov edx,[ang%] cmp edx,1000 jle near loopang inc dword [q%] mov edx,[q%] cmp edx,[dml%] jle near loopq jmp overa .fivhun% DD 500 .overa ret ] endproc rem ============================================================ rem Distribute the wave shape in a circular fashion rem Dummy - never executed def proc_dumcirc for !q%=0 to !dml% !sef%=pmts% + 8*(?t%*!dm% + !q%) |q2%=1.0*!q%/2 for !ang%=0 to 1000 !ex%=int(!cent% + |q2%*cos(!ang%*pi/500)) !ey%=int(!cent% + |q2%*sin(!ang%*pi/500)) !tef%=8*(!ey%*!dm% + !ex%) Fnow(!ex%,!ey%,0)=pnts(!q%,?t%) next next endproc rem ============================================================ rem Assembly Language Routine for Calculating a 3-D distance def proc_norm(opt%) P%=norm% [opt opt% finit fld qword [xd%] fld st0 fmulp st1,st0 fld qword [yd%] fld st0 fmulp st1,st0 fld qword [zd%] fld st0 fmulp st1,st0 faddp st1,st0 faddp st1,st0 fsqrt fstp qword [Ans%] ret ] endproc rem ============================================================ def proc_dumcalc !sef%=pmts% + 8*(?tim%*!dm% + !x%) pnts(!x%,?tim%)=|dk%*!Wht%*|ptts%/!Tot% endproc rem ============================================================ rem Assembly Language Routine for Pnts calculation def proc_calc(opt%) P%=calc% [opt opt% mov eax,0 mov al,[tim%] imul eax,[dm%] add eax,[x%] shl eax,3 mov ebx,pmts% add ebx,eax finit fld qword [dk%] fild dword [Wht%] fmulp st1,st0 fld qword [ptts%] fmulp st1,st0 fild dword [Tot%] fdivp st1,st0 fstp qword [ebx] ret ] endproc rem ============================================================ rem Assembly Language Routine for Master Drawing Routine def proc_Mdrw(opt%) P%=Mdrw% [opt opt% finit fild dword [Wht%] fild dword [ten%] fdivp st1,st0 fstp qword [Sf%] mov edx,[dml%] mov [y%],edx \Initialise y-loop counter .loopyy mov eax,[y%] mov [vy%],eax mov al,4 mov [cd%],al mov eax,[vy%] imul eax,[dm%] shl eax,4 mov [ref%],eax call drw% mov eax,[vy%] dec eax imul eax,[dm%] shl eax,4 call drw% mov al,85 mov [cd%],al call Fdrw% dec dword [y%] mov edx,[y%] cmp edx,1 jge near loopyy jmp overb .ten% DD 10 .overb ret ] endproc rem ============================================================ def proc_dumMdrw |Sf%=!Wht%/10 rem Displays Terrain for !y%=(!dm%-1) to 1 step -1 !vy%=!y% ?cd%=4 !ref%=16*(!vy%*!dm%) call drw% !ref%=16*((!vy%-1)*!dm%) call drw% ?cd%=85 rem Set the colour to a shade of blue according to wave height call Fdrw% next endproc rem ============================================================ rem Set up Arcos array rem 100,000 elements from 0.0 to 1.0 inclusive def proc_arcos_setup local Ans#,lam,a%,pass% dim arcos(100000) for a%=0 to 100000 lam=1.0*(a%/100000) Ans#=1.0*acs(lam) arcos(a%)=Ans# next a% dim rcos% 200, rf% 3 rem Dual-pass assembly, in case of labels for pass%=0 to 2 step 2 proc_rcos(pass%) next pass% cls endproc rem ==================================================================== rem Assembly Language Routine for ArCos Lookup def proc_rcos(opt%) P%=rcos% [opt opt% fld qword [rat%] fild dword [hunth] fmulp st1,st0 fld1 fld1 fld1 faddp st1,st0 fdivp st1,st0 faddp st1,st0 fistp dword [rf%] mov eax,[rf%] shl eax,3 mov ebx,^arcos(0) add ebx,eax fld qword [ebx] fstp qword [Ans%] jmp overt .hunth DD 100000 .overt ret ] endproc rem ==================================================================== rem Assembly Language Routine for ArCos Lookup def proc_alpha(opt%) P%=alpha% [opt opt% mov eax,[y%] imul eax,[dm%] add eax,[x%] shl eax,3 add eax,Fniw% mov [tef%],eax finit fild dword [vht%] \Calc f0r xd fild dword [y%] fild dword [hdis%] faddp st1,st0 fmulp st1,st0 fild dword [hdis%] mov esi,[tef%] fld qword [esi] fild dword [vht%] fsubp st1,st0 fmulp st1,st0 faddp st1,st0 fstp qword [xd%] \Calc f0r xd fild dword [x%] \Calc f0r yd fild dword [dm%] fld1 fld1 faddp st1,st0 fdivp st1,st0 fsubp st1,st0 fld st0 fild dword [vht%] fmulp st1,st0 fchs fstp qword [yd%] \Calc f0r yd fild dword [hdis%] fmulp st1,st0 fchs fstp qword [zd%] \Calc f0r zd call norm% fld qword [xd%] fld qword [Ans%] fdivp st1,st0 fstp qword [rat%] call rcos% mov eax,[y%] imul eax,[dm%] add eax,[x%] shl eax,4 mov [ref%],eax add eax,alpbet% fld qword [Ans%] fstp qword [eax] fild dword [x%] fild dword [dm%] fld1 fld1 faddp st1,st0 fdivp st1,st0 fsubp st1,st0 ftst fstsw ax fstp qword [dmp%] sahf seta al cmp al,1 jne non mov eax,[ref%] add eax,alpbet% fld qword [eax] fchs fstp qword [eax] .non ret ] endproc rem ==================================================================== rem Angle between vertical plane and plane containing position vector def proc_dumalpha !tef%=Fniw% + 8*(!y%*!dm% + !x%) xd=(!vht%*(!y%+!hdis%)+!hdis%*((|(!tef%))-!vht%)) yd=-(!vht%*(!x%-(!dm%/2))) zd=-(!hdis%*(!x%-(!dm%/2))) call norm% |npl%=|Ans% |num%=(!vht%*(!y%+!hdis%)+!hdis%*((|(!tef%))-!vht%)) |rat%=|num%/|npl% call rcos% !ref%=16*(!y%*!dm% + !x%) |(alpbet% + !ref%)=|Ans% if (!x%-(!dm%/2))>0 then |(alpbet% + !ref%)=-|(alpbet% + !ref%) endproc rem----------------------------------------------------------------------- rem Assembly Language Routine for ArCos Lookup def proc_beta(opt%) P%=beta% [opt opt% finit fild dword [x%] fild dword [dm%] fld1 fld1 faddp st1,st0 fdivp st1,st0 fsubp st1,st0 fstp qword [xd%] fild dword [y%] fild dword [hdis%] faddp st1,st0 fstp qword [yd%] mov esi,[tef%] fld qword [esi] fild dword [vht%] fsubp st1,st0 fstp qword [zd%] call norm% fild dword [y%] fild dword [hdis%] fadd st1,st0 fmulp st1,st0 mov esi,[tef%] fld qword [esi] fild dword [vht%] fsubr st1,st0 fmulp st1,st0 faddp st1,st0 fld qword [Ans%] fdivp st1,st0 fld qword [Rr%] fdivp st1,st0 fstp qword [rat%] call rcos% mov eax,[ref%] add eax,alpbet% fld qword [eax] fld st0 fsin fchs add eax,8 fld qword [Ans%] fst qword [eax] fmulp st1,st0 fstp qword [AX%] fcos fld qword [Ans%] fmulp st1,st0 fstp qword [AY%] mov eax,[ref%] mov ebx,prs% add ebx,eax fld qword [AX%] fstp qword [ebx] fld qword [AY%] add ebx,8 fstp qword [ebx] ret ] endproc rem ==================================================================== rem 2-D Projection View-Angle def proc_dumbeta |xd%=(!x%-(!dm%/2)) : |yd%=(!y%+!hdis%) : |zd%=((|(!tef%))-!vht%) call norm% rdis=|Ans% yfact=(!hdis%*(!hdis%+!y%)) : zfact=(!vht%*(!vht%-(|(!tef%)))) |rat%=(yfact+zfact)/(rdis*|Rr%) call rcos% !ref%=16*(!y%*!dm% + !x%) |(alpbet% + !ref% + 8)=|Ans% rem Convert from Polar to x/y Coords for Printing to Screen |AX%=|(alpbet% + !ref% + 8)*-sin(|(alpbet% + !ref%)) |AY%=|(alpbet% + !ref% + 8)*cos(|(alpbet% + !ref%)) rem Place AX and AY at consecutive places @ !ref% in the array prs% call Floc% endproc rem =====================================================================