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 =====================================================================