rem Field Assembly ..... Rev 8.0
      rem A J Tooth / 25th July 2003
      rem Modified with Assembly Routines April 2004
      rem Default background picture incorporated in compiled version // July 2005
      on error if (err=17) then quit
      
      rem Call Utility Procedures
      proc_fullscreen  :  rem Set up use of Full Screen
      himem=lomem + 200000000
      *FLOAT64
      
      rem Setup parameters
      proc_setup
      
      rem Randomize landscape
      mouse on 2 : colour 3
      print tab(15,15);"STATUS: Randomising Landscape       "
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &100
      proc_land(1,1)
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &20
      
      rem Smooth the landscape
      ps%=0
      while ps%<Ps%
        print tab(15,15);"STATUS: Smoothing Pass ";ps%+1;"           "
        sys "GetCurrentProcess" to hprocess%
        sys "SetPriorityClass", hprocess%, &100
        proc_smooth
        sys "GetCurrentProcess" to hprocess%
        sys "SetPriorityClass", hprocess%, &20
        
        ps%+=1
      endwhile
      rem Resets all areas of WATER to level -1
      proc_water
      
      rem Project landscape plane onto viewing screen
      print tab(15,15);"STATUS: Calculating Perspective       "
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &100
      proc_pers
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &20
      
      rem Draw the landscape
      mouse off : cls : clg
      proc_Fdraw(0)
      
      rem Press any key to proceed to light section
      sp$=get$
      mouse on 2
      
      rem Picks a random direction for the sunlight to come from
      proc_lightran
      
      print tab(15,0);"STATUS: Calculating Normals       "
      
      rem Calculates the normals to each of the surface elements
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &100
      proc_norm
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &20
      
      clg : mouse off
  680 clg : proc_Fdraw(1)
      
      rem Press any key to exit
      sp$=get$
      if sp$="h" then !vt%=int(1.5*!vt%) : proc_pers : goto 680
      if sp$="l" then !vt%=int(0.5*!vt%) : proc_pers : goto 680
      if sp$="f" then !hdis%=int(0.75*!hdis%) : proc_pers : goto 680
      if sp$="b" then !hdis%=int(1.25*!hdis%) : proc_pers : goto 680
      if sp$=" " then run
      quit
      end
      rem End of Programme --------------------------------------------------
      
      rem Setup parameters
      def proc_setup
      dim vt% 3, hdis% 3, dm% 3, Siz% 3, scl% 3
      
      colour 3
      *FONT Verdana,12,B
      
      rem Place Origin and Draw Axes
      input " Vertical Height (400)",!vt%
      if !vt%=0 then !vt%=400
      input " Horizontal Distance (1500)",!hdis%
      if !hdis%= 0 then !hdis%=1500
      input " Ultra-High, High or Low Resolution? ( -u- or -h- or -l-)(l)",in$
      if in$="" then in$="l"
      if in$="h" then !vt%=2*!vt% : !hdis%=2*!hdis%
      input " How many smoothing passes (7)",Ps%
      if Ps%=0 then Ps%=7
      
      print:print:print " Left-click to use default jpeg for background."
      print " Otherwise right-click to select an alternative picture."
      
      rem Ensure debounce taken care of
      repeat
        mouse x,y,a&
      until a&=0
      
      repeat
        mouse x,y,a&
      until a&=4 or a&=1
      
      if a&=4 then
        path$=@dir$
        Pname$="Sky3.jpg"
      else
        proc_locload
      endif
      
      rem Set up square field of side dm% (256, 512 or 1024)
      case in$ of
        when "u" : !dm%=1024 : !scl%=7000
        when "h" : !dm%=512  : !scl%=15000
        when "l" : !dm%=256  : !scl%=15000
      endcase
      
      !Siz%=(!dm%+1)*(!dm%+1)
      
      rem Dimension Field Array and Projection Array
      print tab(15,15);"STATUS: Dimensioning Arrays       "
      
      dim alpbet(!dm%,!dm%,2)
      dim ps(!dm%,!dm%,2)
      dim trac%(75000,6)
      dim Norm(!dm%,!dm%,3)
      dim work% 3, r% 0,g% 0,b% 0,Red% 0,Gre% 0,Blu% 0
      dim xmax% 3, zmax% 3, fset% 3, R% 7
      dim xl% 3,yl% 3,xh% 3,yh% 3,vy% 3,par1% 7,par2% 7,par3% 7,par4% 7
      dim psx% 3, psy% 3,dmp% 3,x% 3,y% 3,med% 7,yfix% 3,Ref% 3, iRef1% 3, iRef2% 3
      dim lvx% 3, lvy% 3, nine% 7, xtem% 3, ytem% 3, lx% 0, ly% 0
      dim u1% 7, u2% 7, u3% 7, m1% 7, k% 3, dum% 0
      dim xd% 7, yd% 7, zd% 7, num% 7, npl% 7, yfact% 7, zfact% 7, rdis% 7
      dim frstx% 3, frsty% 3, cntx% 3, cnty% 3, scndx% 3, scndy% 3
      dim xmid% 3, ymid% 3, xnew% 3, ynew% 3, avg% 7, pert% 7, fct% 7, fact% 7, red% 7
      dim Svec1% 7, Svec2% 7, Svec3% 7, Tes1% 3, Tes2% 3, Tes3% 7, TesH% 7
      dim Term1% 7, Term2% 7, Term3% 7, ResCos% 7, fac% 7
      
      dim Fnow% 8*!Siz%, Fnowt% !Siz%, drw% 1000, Fsm% 8*!Siz%, xfer% 1000, smth% 1000
      dim alph% 1000, Ref8% 3, p13% 7, p23% 7, p33% 7, beth% 1000
      
      colour 7
      rem Centres Graphics and sets scale
      !xmax%=1024 : !zmax%=768
      !fset%=0
      |fct%=0.5*rnd(1) : rem Used in Perturbation Procedure
      ?dum%=40+rnd(160)
      
      rem Dual-pass assembly, in case of labels
      for i%=0 to 2 step 2
        proc_drw(i%)
        proc_xfer(i%)
        proc_smth(i%)
        proc_alph(i%)
        proc_beth(i%)
      next i%
      
      origin !xmax%,!fset%
      |R%=sqr(!vt%*!vt%+!hdis%*!hdis%)
      endproc
      rem----------------------------------------------------------------------
      
      rem Locate and load picture into memory
      def proc_locload
      local g%
      
      dim pf% 75, ff% 18, fm% 255
      !pf% = 76
      pf%!4 = @hwnd%
      pf%!12 = ff%
      pf%!28 = fm%
      pf%!32 = 256
      pf%!52 = 6
      $ff% ="JPG files"+chr$0+"*.jpg"+chr$0+chr$0
      sys "GetOpenFileName", pf% to result%
      if result%<>0 then fullname$ = fn_nulterm$(fm%)
      
      rn%=len(fullname$)
      g%=0 : file$=""
      repeat
        n$=mid$(fullname$,rn%-g%,1)
        if n$<>"\" then file$=n$+file$ : g%+=1
      until n$="\"
      Pname$=file$
      
      path$=left$(fullname$,rn%-g%)
      rem Change the current Directory
      command$="CD "+chr$(34)+path$+chr$(34)
      oscli command$
      endproc
      rem ===============================================================================================
      
      rem Convert null-terminated string to BASIC format
      def fn_nulterm$(Z%)
      local A$
      while ?Z%<>0
        A$+=chr$(?Z%)
        Z%+=1
      endwhile
      =A$
      rem ===============================================================================================
      
      rem Draw Points on Screen
      def proc_Fdraw(Flg%)
      local a$,xhorL%,xhorR%,diff%,g%
      
      rem Set up SKY as Backdrop
      pic$=path$+Pname$
      
      edge%=int(!scl%*ps(0,!dm%,2))
      rem Printss the above JPG file to the screen.
      proc_pic(edge%,pic$)
      
      rem Displays Terrain
      cld%=0
      for !y%=!dm% to 1 step -1
        !vy%=!y%
        move !scl%*ps(0,!vy%,1),!scl%*ps(0,!vy%,2)
        move !scl%*ps(0,!vy%-1,1),!scl%*ps(0,!vy%-1,2)
        for !x%=1 to !dm%
          Ref8%=8*!x% + 8*(!dm%+1)*!vy%
          indRef8%=Fnow% + Ref8%
          ch%=int(|indRef8%)
          
          if ch%>50 then sn%=1 else sn%=0
          if ch%>0 then
            col%=100-2*ch% : if col%<0 then col%=0
          else
            col%=ch%
          endif
          rem If col%<0 then set flag to indicate below water level
          rem Below water level, colour is naturally BLUE
          if col%<0 then t%=1 else t%=0
          rem Otherwise, colour ranges from green to dark brown with altitude
          rem On second pass, colours vary according to the light intensity
          if Flg%=1 then
            rem If it's land, shade it. If it's water, reflect the background in it.
            case t% of
              when 0 : |fac%=0.00001 + 0.5*fn_cosPh
              when 1 : |fac%=0.00001 : rgb%=fn_Refl
            endcase
          else
            |fac%=0.00001
          endif
          
          rem Set r/g/b levels
          if (Flg%=1 and t%=1) then
            if rgb%=-1 then rgb%=0
            !work%=rgb%
            ?Red%=!work% mod 256 : !work%=!work%-?Red% : !work%=!work%/256
            ?Gre%=!work% mod 256 : !work%=!work%-?Gre% : !work%=!work%/256
            ?Blu%=!work% mod 256
            f=0.9
            ?Red%=int(f*?Red%) : ?Gre%=int(f*?Gre%) : ?Blu%=int(f*?Blu%)
          else
            rem If sn%=1 then there's snow, so choose WHITE
            if sn%=0 then
              ?r%=30 : ?g%=5+col% : ?b%=30
            else
              ?r%=127 : ?g%=127 : ?b%=127
            endif
            ?Red%=int((?dum%*!vy% + ?r%*(1.3+|fac%)*(1-t%)*(!dm%-!vy%))/!dm%)
            ?Gre%=int((?dum%*!vy% + ?g%*(1.3+|fac%)*(1-t%)*(!dm%-!vy%))/!dm%)
            ?Blu%=int((?dum%*!vy% + ?b%*(1.3+|fac%)*(!dm%-!vy%))/!dm%)
          endif
          
          |par1%=ps(!x%,!vy%,1)
          |par2%=ps(!x%,!vy%,2)
          !vy%-=1
          |par3%=ps(!x%,!vy%,1)
          |par4%=ps(!x%,!vy%,2)
          !vy%+=1
          rem Call Assembly Routine to change colours and plot triangles
          call drw%
          
        next
      next
      endproc
      
      rem --------------------------------------------------
      
      rem FUNCTIONS
      rem =========
      
      rem Calculates Cosine of Light Incidence to get Intensity
      def fn_cosPh
      local hit%
      
      Ref8%=8*!x% + 8*(!dm%+1)*!vy%
      indRef8%=Fnow% + Ref8%
      
      rem Check if the sunlight hits an intervening hill
      !k%=0 : hit%=0
      repeat
        !k%+=25
        rem Calculate test vector point along line of sunlight
        !Tes1%=int(!x% - !k%*|Svec1%)
        !Tes2%=int(!vy% - !k%*|Svec2%)
        |Tes3%=(|indRef8% - !k%*|Svec3%)
        
        if (!Tes1%>-1 and !Tes1%<(!dm%+1) and !Tes2%>-1 and !Tes2%<(!dm%+1)) then
          Ref8%=8*!Tes1% + 8*(!dm%+1)*!Tes2%
          indRef8%=Fnow% + Ref8%
          |TesH%=|indRef8%
          if |Tes3%<|TesH% then hit%=1
        endif
        
      until (|Tes3%>=300.0 or hit%=1)
      
      |Term1%=|Svec1%*Norm(!x%-1,!vy%,1)
      |Term2%=|Svec2%*Norm(!x%-1,!vy%,2)
      |Term3%=|Svec3%*Norm(!x%-1,!vy%,3)
      
      rem ResCos is the resultant Cosine of the angle of incidence
      |ResCos%=|Term1% + |Term2% + |Term3%
      rem If the sunlight hit a hill, colour to be
      rem no brighter than for 90Deg incidence angle.
      if hit%=1 then if |ResCos%>0.00001 then |ResCos%=0.00001
      =|ResCos%
      
      rem ----------------------------------------------------------
      
      rem PROCEDURES
      rem ==========
      rem Calculates perspective
      def proc_pers
      
      print tab(15,15);"STATUS: Calculating Perspective       "
      |R%=sqr(!vt%*!vt%+!hdis%*!hdis%)
      for !x%=0 to !dm%
        for !y%=0 to !dm%
          skp%=0
          on error local proc_error
          if skp%=1 then goto 3420
          
          proc_alpha
          proc_beta
          
 3420   next
      next
      endproc
      
      rem-----------------------------------------------------
      
      rem Angle between vertical plane and plane containing position vector
      def proc_alpha
      
      rem Various mundane calculations
      call alph%
      
      if |npl%=0 then alpbet(!x%,!y%,1)=pi/2 else alpbet(!x%,!y%,1)=acs(|num%/|npl%)
      
      if (!x%-(!dm%/2))>0 then alpbet(!x%,!y%,1)=-alpbet(!x%,!y%,1)
      endproc
      
      rem-----------------------------------------------------
      
      rem 2-D Projection View-Angle
      def proc_beta
      
      call beth%
      
      alpbet(!x%,!y%,2)=acs((|yfact%+|zfact%)/(|rdis%*|R%))
      
      rem Convert from Polar to x/y Coords for Printing to Screen
      ps(!x%,!y%,1)=alpbet(!x%,!y%,2)*-sin(alpbet(!x%,!y%,1))
      ps(!x%,!y%,2)=alpbet(!x%,!y%,2)*cos(alpbet(!x%,!y%,1))
      endproc
      
      rem --------------------------------------------------
      
      rem Utility Procedures
      rem ==================
      
      rem Set up use of Full Screen
      def proc_fullscreen
      sys "GetSystemMetrics", 0 to xscreen%
      sys "GetSystemMetrics", 1 to yscreen%
      sys "SetWindowLong",@hwnd%,-16,&16000000
      sys "SetWindowPos",@hwnd%,-1,0,0,xscreen%,yscreen%,0
      vdu 23,22,xscreen%;yscreen%;8,16,16,1 : rem Set fullscreen mode
      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 Randomizes landscape
      def proc_land(Lps%,Flg%)
      local nd%,ind%,Max%
      
      Max%=350
      stp%=0
      nd%=0
      
      rem Flg%=1 indicates the 1st pass
      if Flg%=1 then
        rem Set end-points
        !xl%=0 : !yl%=0 : !xh%=!dm% : !yh%=!dm%
        !Ref%=!xl% + (!dm%+1)*!yl%
        Fnowt%?!Ref%=1
        !Ref%=!xh% + (!dm%+1)*!yh%
        Fnowt%?!Ref%=1
        !Ref%=!xl% + (!dm%+1)*!yh%
        Fnowt%?!Ref%=1
        !Ref%=!xh% + (!dm%+1)*!yl%
        Fnowt%?!Ref%=1
      endif
      
      nd%=0 : ind%=0
      lp%=Lps% : !frstx%=!xl% : !frsty%=!yl%
      
      repeat
        !cntx%=!frstx%+1 : !cnty%=!frsty%+1
        !scndx%=-1 : !scndy%=-1
        rem Find the next perturbed point
        repeat
          !Ref%=!cntx% + (!dm%+1)*!frsty%
          if Fnowt%?!Ref%=1 then !scndx%=!cntx% else !cntx%+=1
        until (!cntx%=(!dm%+1) or !scndx%>0)
        repeat
          !Ref%=!frstx% + (!dm%+1)*!cnty%
          if Fnowt%?!Ref%=1 then !scndy%=!cnty% else !cnty%+=1
        until (!cnty%=(!dm%+1) or !scndy%>0)
        
        rem Set exit criterion
        if (((!scndx%-!frstx%)>1) and ((!scndy%-!frsty%)>1)) then
          proc_midl(Flg%)
          if (!scndx%-!frstx%>2) then
            ind%+=1
            trac%(ind%,1)=!xnew%  : trac%(ind%,2)=!ynew%
            trac%(ind%,3)=!frstx% : trac%(ind%,4)=!frsty%
            trac%(ind%,5)=!scndx% : trac%(ind%,6)=!scndy%
          else
            ind%=0
          endif
        endif
        
        if (!frstx%=(!dm%-2) and !scndx%=!dm% and !frsty%=(!dm%-2) and !scndy%=!dm%) then nd%=1 : goto 4720
        rem Move reference points along
        if !scndx%=!dm% then
          !frstx%=0
          if !scndy%=!dm% then
            !frsty%=0 : lp%+=1
          else
            !frsty%=!scndy%
          endif
        else
          !frstx%=!scndx%
        endif
        
        rem Create new corner points
        if (!frstx%=0 and !frsty%=0 and ind%>0) then
          for f%=1 to ind%
            proc_fill(f%)
          next f%
          ind%=0
        endif
        Flg%=2
 4720 until (nd%=1)
      endproc
      
      rem +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      rem Set water values
      def proc_water
      for !x%=0 to !dm%
        for !y%=0 to !dm%
          Ref8%=8*!x% + 8*(!dm%+1)*!y%
          indRef8%=Fnow% + Ref8%
          if |indRef8%<0.0 then |indRef8%=-1.0
        next
        Ref8%=8*!x% + 8*(!dm%+1)*!dm%
        indRef8%=Fnow% + Ref8%
        |indRef8%=-1.0
      next
      endproc
      
      rem +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      rem Perturbs mid-point
      def proc_midl(Fl%)
      
      rem Reduce max perturbation with each loop.
      |fact%=1.1+|fct%
      |red%=(|fact%)^(1.2*lp%)
      !xmid%=int(!frstx%+((!scndx%-!frstx%)/2)) : !ymid%=int(!frsty%+((!scndy%-!frsty%)/2))
      
      rem If there ARE still mid-points, select a mid-point
      if Fl%=1 then
        sw%=1
      else
        sw%=(2*(rnd(2)-1))-1
      endif
      |pert%=0.0 + 1.00001*sw%*rnd(int(Max%/|red%))
      !Ref%=!xmid% + (!dm%+1)*!ymid%
      Fnowt%?!Ref%=1
      
      |avg%=0.0
      Ref8%=8*!frstx% + 8*(!dm%+1)*!frsty%
      indRef8%=Fnow% + Ref8%
      |avg%+=|indRef8%
      Ref8%=8*!frstx% + 8*(!dm%+1)*!scndy%
      indRef8%=Fnow% + Ref8%
      |avg%+=|indRef8%
      Ref8%=8*!scndx% + 8*(!dm%+1)*!frsty%
      indRef8%=Fnow% + Ref8%
      |avg%+=|indRef8%
      Ref8%=8*!scndx% + 8*(!dm%+1)*!scndy%
      indRef8%=Fnow% + Ref8%
      |avg%+=|indRef8%
      
      |avg%=(|avg%)/4
      
      Ref8%=8*!xmid% + 8*(!dm%+1)*!ymid%
      indRef8%=Fnow% + Ref8%
      
      |indRef8%=|avg% + |pert%
      
      !xnew%=!xmid% : !ynew%=!ymid%
      endproc
      
      rem ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      rem Fill in new points
      def proc_fill(Ind%)
      local xf%,yf%,xn%,yn%,xs%,ys%
      
      xn%=trac%(Ind%,1) : yn%=trac%(Ind%,2)
      xf%=trac%(Ind%,3) : yf%=trac%(Ind%,4)
      xs%=trac%(Ind%,5) : ys%=trac%(Ind%,6)
      
      rem Set heights of new points
      rem Uses quasi-cubic interpolation
      Ref8%=8*(xn%-1) + 8*(!dm%+1)*(yn%-1)
      indRef8%=Fnow% + Ref8%
      
      !Ref%=xn% + (!dm%+1)*yf%
      Ref8%=8*xn% + 8*(!dm%+1)*yf%
      indRef8%=Fnow% + Ref8%
      Ref8a%=8*xf% + 8*(!dm%+1)*yf%
      indRef8a%=Fnow% + Ref8a%
      Ref8b%=8*xs% + 8*(!dm%+1)*yf%
      indRef8b%=Fnow% + Ref8b%
      Fnowt%?!Ref%=1 : |indRef8%=(|indRef8a% + |indRef8b%)/2
      
      !Ref%=xf% + (!dm%+1)*yn%
      Ref8%=8*xf% + 8*(!dm%+1)*yn%
      indRef8%=Fnow% + Ref8%
      Ref8a%=8*xf% + 8*(!dm%+1)*yf%
      indRef8a%=Fnow% + Ref8a%
      Ref8b%=8*xs% + 8*(!dm%+1)*ys%
      indRef8b%=Fnow% + Ref8b%
      Fnowt%?!Ref%=1 : |indRef8%=(|indRef8a% + |indRef8b%)/2
      
      !Ref%=xn% + (!dm%+1)*ys%
      Ref8%=8*xn% + 8*(!dm%+1)*ys%
      indRef8%=Fnow% + Ref8%
      Ref8a%=8*xf% + 8*(!dm%+1)*ys%
      indRef8a%=Fnow% + Ref8a%
      Ref8b%=8*xs% + 8*(!dm%+1)*ys%
      indRef8b%=Fnow% + Ref8b%
      Fnowt%?!Ref%=1 : |indRef8%=(|indRef8a% + |indRef8b%)/2
      
      !Ref%=xs% + (!dm%+1)*yn%
      Ref8%=8*xs% + 8*(!dm%+1)*yn%
      indRef8%=Fnow% + Ref8%
      Ref8a%=8*xs% + 8*(!dm%+1)*yf%
      indRef8a%=Fnow% + Ref8a%
      Ref8b%=8*xs% + 8*(!dm%+1)*ys%
      indRef8b%=Fnow% + Ref8b%
      Fnowt%?!Ref%=1 : |indRef8%=(|indRef8a% + |indRef8b%)/2
      
      endproc
      
      rem --------------------------------------------------------------
      
      rem Smooths the Field
      def proc_smooth
      
      !yfix%=8*(!dm%+1) : |nine%=9.0
      
      !x%=0 : Ref1%=0 : Ref2%=8*!dm% : Ref3%=0 : Ref4%=!yfix%*!dm%
      repeat
        indRef1%=Fsm% + Ref1%
        |indRef1%=0
        Ref1%+=!yfix%
        
        indRef2%=Fsm% + Ref2%
        |indRef2%=0
        Ref2%+=!yfix%
        
        indRef3%=Fsm% + Ref3%
        |indRef3%=0
        Ref3%+=8
        
        indRef4%=Fsm% + Ref4%
        |indRef4%=0
        Ref4%+=8
        
        !x%+=1
      until (!x%=(!dm%+1))
      
      for !y%=1 to (!dm%-1)
        for !x%=1 to (!dm%-1)
          
          rem Call Assembly Routine for Smoothing
          call smth%
          
        next
      next
      
      rem Calls Assembly Routine for Transfer
      call xfer%
      
      endproc
      
      rem +++++++++++++++++++++++++++++++++++++++++++++++++++
      
      rem Random Light Source
      def proc_lightran
      local Ph,Th
      
      Ph=rnd(45) : Th=rnd(90)-45
      |Svec1%=-(cos(rad(Ph)))*cos(rad(Th))
      |Svec2%=-(cos(rad(Ph)))*sin(rad(Th))
      |Svec3%=-sin(rad(Ph))
      
      endproc
      
      rem +++++++++++++++++++++++++++++++++++++++++++++++++++
      
      rem Calculate Normals to every triangulation
      def proc_norm
      
      for !y%=!dm% to 1 step -1
        for !x%=0 to (!dm%-1)
          
          Ref8%=8*!x% + 8*(!dm%+1)*!y%
          indRef8%=Fnow% + Ref8%
          
          |p13%=|indRef8%
          
          indRef8%=Fnow% + Ref8% - 8*(!dm%+1)
          |p23%=|indRef8%
          
          indRef8%=Fnow% + Ref8% + 8
          |p33%=|indRef8%
          
          |u1%=|p13%-|p33% : |u2%=|p23%-|p13%
          
          |m1%=sqr(|u1%*|u1% + |u2%*|u2% + 1.0)
          
          Norm(!x%,!y%,1)=-|u1%/|m1%
          Norm(!x%,!y%,2)=-|u2%/|m1%
          Norm(!x%,!y%,3)=-|u3%/|m1%
          
        next
      next
      endproc
      
      rem +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      def proc_pic(Edge%,Pic$)
      local nel%,ps%
      
      ps%=0
      picture$ = Pic$ : rem. Image filename
      
 6810 xpos% = 0                  : rem. X position (pixels)
      ypos% = 0                  : rem. Y position (pixels)
      xsize% = 1024              : rem. Width (pixels)
      ysize% = 768-(Edge%/2)     : rem. Height (pixels)
      
      sys "LoadLibrary", "OLEAUT32.DLL"
      sys "GetModuleHandle", "OLEAUT32.DLL" to oleaut32%
      sys "GetProcAddress", oleaut32%, "OleLoadPicturePath" to olpp%
      if olpp%=0 error 0, "Could not get address of OleLoadPicturePath"
      
      dim iid% 15, gpp% 3, hmw% 3, hmh% 3, picture% 513
      sys "MultiByteToWideChar", 0, 0, picture$, len(picture$), picture%, 256 to nel%
      picture%!(2*nel%) = 0
      
      iid%!0  = &7BF80980
      iid%!4  = &101ABF32
      iid%!8  = &AA00BB8B
      iid%!12 = &AB0C3000
      
      sys olpp%, picture%, 0, 0, 0, iid%, gpp%
      ps%+=1
      if (ps%=1 and !gpp% = 0) then
        lg%=len(picture$)
        picture$=left$(picture$,(lg%-4)) + "JPG"
      endif
      if (ps%=1 and !gpp%=0) goto 6810
      
      if (ps%=2 and !gpp%=0) then quit
      
      sys !(!!gpp%+24), !gpp%, hmw% : rem. IPicture::get_Width
      sys !(!!gpp%+28), !gpp%, hmh% : rem. IPicture::get_Height
      sys !(!!gpp%+32), !gpp%, @memhdc%, xpos%, ypos%, xsize%, ysize%, 0, !hmh%, !hmw%, -!hmh%, 0 to res%
      if res% error 0, "IPicture::Render failed"
      
      sys !(!!gpp%+8), !gpp% : rem. IPicture::Release
      sys "InvalidateRect", @hwnd%, 0, 0
      sys "UpdateWindow", @hwnd%
      
      endproc
      
      rem ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      rem Reflection Routine
      def fn_Refl
      local Rgb%,lam,hit%,Mq,q1,q2,q3,Dx%
      
      Dx%=!x%-(!dm%/2)
      
      rem (q1,q2,q3) is the line-of-sight vector to the point on the SEA ONLY.
      q1=Dx% : q2=!vy%+!hdis% : q3=-!vt%
      Mq=sqr(q1*q1 + q2*q2 + q3*q3)
      rem (s1,s2,s3) is the UNIT reflected vector FROM the point on the sea.
      s1=Dx%/Mq : s2=(!vy%+!hdis%)/Mq : s3=!vt%/Mq
      
      rem lam is the factor determining where the reflection vector meets the sky
      rem com is a common expression cropping up in the calculations.
      com=(!dm%-!vy%)/(!vy%+!hdis%)
      lam=(!dm%-!vy%)*Mq/(!vy%+!hdis%)
      
      rem Does the reflection hit intervening land, or carry on to meet the sky?
      !k%=0
      repeat
        !k%+=25
        !Tes1%=!x%+int(!k%*s1) : !Tes2%=!vy%+int(!k%*s2) :|Tes3%=!k%*s3
        if (!Tes1%>-1 and !Tes1%<(!dm%+1) and !Tes2%>-1 and !Tes2%<(!dm%+1)) then
          Ref8%=8*!Tes1% + 8*(!dm%+1)*!Tes2%
          indRef8%=Fnow% + Ref8%
          if |indRef8%>|Tes3% then hit%=1
        endif
        if !k%>lam then hit%=2
        rem This routine leaves with the appropriate value of k%
      until (hit%=1 or hit%=2)
      
      rem Work out relevant angles
      proc_alp(!Tes1%,!Tes2%,com)
      proc_bet(!Tes1%,!Tes2%,com)
      rem Find out what colour the sky is at the point where
      rem the reflection vector "hits" the sky on screen.
      Rgb%=tint(posx,posy)
      =Rgb%
      
      rem-----------------------------------------------------
      
      rem Modified Functional Procedures specific to Reflection Routine
      rem =============================================================
      rem Angle between vertical plane and plane containing position vector
      def proc_alp(Xx%,Yy%,Com)
      local oda,num,xd,yd,zd
      
      rem Calculates modulus of vector from origin to
      rem where the reflection "hits" the sky.
      xd=(-!hdis%*!vt%*Com) - (!dm%*!vt%)
      yd=Dx%*!vt%*(1+Com)
      zd=Dx%*!hdis%*(1+Com)
      oda=sqr(xd*xd+yd*yd+zd*zd)
      
      rem Works out the angle between the plane containing the line of sight
      rem vector and the vertical plane.
      num=(!vt%*!hdis%*Com) + (!dm%*!vt%)
      if oda=0 then alp=pi/2 else alp=acs(num/oda)
      
      if (Xx%-(!dm%/2))>0 then alp=-alp
      endproc
      
      rem-----------------------------------------------------
      
      rem 2-D Projection View-Angle
      def proc_bet(Xx%,Yy%,Com)
      local rdis,xd,yd,zd,yfact,zfact
      
      rem Calculates modulus of line-of-sight vector to
      rem where the reflection vector hits the sky.
      xd=Dx%*(1+Com)
      yd=!dm%+!hdis%
      zd=!vt%*(Com-1)
      rdis=sqr(xd*xd+yd*yd+zd*zd)
      
      rem Works out the angle between the line-of-sight vector to the origin
      rem and the line-of-sight vector to where the reflection hits the sky
      yfact=!hdis%*(!hdis%+!dm%) : zfact=!vt%*!vt%*(1-Com)
      bet=acs((yfact+zfact)/(rdis*|R%))
      
      rem Convert from Polar to x/y Co-ords to identify where on the screen
      rem to fetch the RGB value of the screen at that point.
      posx=!scl%*bet*-sin(alp)
      posy=!scl%*bet*cos(alp)
      endproc
      rem----------------------------------------------------------------------
      
      rem Assembly Routine for Plotting
      def proc_drw(opt%)
      P%=drw%
      
      [opt opt%
      mov al,[vdsetcol%]      ;Changes the rgb vlue f0r colr 10
      call "oswrch"
      mov al,[cl%]
      call "oswrch"
      mov al,[paltyp%]
      call "oswrch"
      mov al,[Red%]
      call "oswrch"
      mov al,[Gre%]
      call "oswrch"
      mov al,[Blu%]
      call "oswrch"
      
      mov al,[vdgcol%]        ;Calls gc0l routine
      call "oswrch"
      mov al,[md%]
      call "oswrch"
      mov al,[cl%]
      call "oswrch"
      
      finit
      fild dword [scl%]
      fld qword [par1%]
      fmul st0,st1
      fistp dword [psx%]
      fld qword [par2%]
      fmul st0,st1
      fistp dword [psy%]
      fistp dword [dmp%]
      
      mov al,[pltcode%]      ;Calls pl0t 85 routine
      call "oswrch"
      mov al,[tri%]
      call "oswrch"
      mov bx,[psx%]
      mov al,bl
      call "oswrch"
      mov al,bh
      call "oswrch"
      mov bx,[psy%]
      mov al,bl
      call "oswrch"
      mov al,bh
      call "oswrch"
      
      fild dword [scl%]
      fld qword [par3%]
      fmul st0,st1
      fistp dword [psx%]
      fld qword [par4%]
      fmul st0,st1
      fistp dword [psy%]
      fistp dword [dmp%]
      
      mov al,[pltcode%]      ;Calls pl0t 85 routine
      call "oswrch"
      mov al,[tri%]
      call "oswrch"
      mov bx,[psx%]
      mov al,bl
      call "oswrch"
      mov al,bh
      call "oswrch"
      mov bx,[psy%]
      mov al,bl
      call "oswrch"
      mov al,bh
      call "oswrch"
      
      jmp round              ;Various constants
      .vdsetcol% db 19
      .paltyp% db 16
      .vdgcol% db 18
      .md% db 0
      .cl% db 10
      .pltcode% db 25
      .tri% db 85
      
      .round ret
      ]
      endproc
      
      rem-----------------------------------------------------
      
      rem Assembly Routine for Number Transfer
      def proc_xfer(opt%)
      P%=xfer%
      
      [opt opt%
      
      mov eax,0             ;Initialise loop variables
      mov [Ref%],eax
      mov [lvx%],eax
      mov [lvy%],eax
      
      .loop
      mov eax,[Ref%]
      mov ebx,Fsm%[eax]     ;Shift Fsm t0 Fnow
      mov Fnow%[eax],ebx
      add eax,4
      mov ebx,Fsm%[eax]     ;Needs 2 attempts because the numbers are 64-bit
      mov Fnow%[eax],ebx
      add eax,4
      mov [Ref%],eax
      inc dword [lvx%]
      mov eax,[lvx%]
      cmp eax,[dm%]
      jle loop              ;Repeat If lvx%<=dm%
      
      mov eax,0             ;Revert lvx% t0 0
      mov [lvx%],eax
      inc dword [lvy%]
      mov eax,[lvy%]
      cmp eax,[dm%]
      jle loop              ;Repeat If lvy%<=dm%
      
      ret
      ]
      endproc
      
      rem Assembly Routine for Smoothing
      def proc_smth(opt%)
      P%=smth%
      
      [opt opt%
      
      finit
      fldz
      fstp qword [med%]    ;Initialise med% t0 zero
      mov al,0
      mov [lx%],al
      mov [ly%],al
      
      .rep
      mov eax,[x%]
      dec eax
      mov [xtem%],eax
      mov eax,0
      mov al,[lx%]
      add [xtem%],eax
      mov eax,[xtem%]
      shl eax,3
      mov [Ref%],eax
      mov eax,[y%]
      dec eax
      mov [ytem%],eax
      mov eax,0
      mov al,[ly%]
      add [ytem%],eax
      mov eax,[ytem%]
      mov ebx,[yfix%]
      imul eax,ebx
      add [Ref%],eax        ;Ref%= 8*(x%-1+lx%) + yfix%*(y%-1+ly%) at this point
      
      mov eax,[Ref%]
      fld qword [med%]
      fld qword Fnow%[eax]
      faddp st1,st0
      fstp qword [med%]
      
      inc byte [lx%]
      mov al,[lx%]
      cmp al,2
      jle near rep
      
      mov al,0
      mov [lx%],al
      inc byte [ly%]
      mov al,[ly%]
      cmp al,2
      jle near rep
      
      mov eax,[x%]
      shl eax,3
      mov [Ref%],eax
      mov eax,[y%]
      mov ebx,[yfix%]
      imul eax,ebx
      add [Ref%],eax
      mov eax,[Ref%]
      
      fld qword [med%]
      fld qword [nine%]
      fdivp st1,st0
      fstp qword Fsm%[eax]
      
      ret
      ]
      endproc
      
      rem Assembly Routine for Alpha Calc.
      def proc_alph(opt%)
      P%=alph%
      
      [opt opt%
      
      mov eax,[dm%]
      inc eax
      mov ebx,[y%]
      imul eax,ebx
      mov ebx,[x%]
      add eax,ebx
      shl eax,3
      mov [Ref8%],eax       ;Calculates Ref8%, the Offset in the array Fnow%
      
      finit
      fld qword Fnow%[eax]
      fild dword [vt%]
      fsubp st1,st0
      fild dword [hdis%]
      fmul st1,st0
      fild dword [y%]
      faddp st1,st0
      fild dword [vt%]
      fmulp st1,st0
      faddp st1,st0
      fstp qword [xd%]      ;Calculates xd%
      
      fild dword [dm%]
      fld1
      fld1
      faddp st1,st0
      fdivp st1,st0
      fild dword [x%]
      fsubp st1,st0         ;dm%/2 - x% calculated here
      fld st0               ;... And is duplicated
      
      fild dword [vt%]
      fmulp st1,st0
      fstp qword [yd%]      ;Calculates yd%
      fild dword [hdis%]
      fmulp st1,st0
      fst qword [zd%]       ;Calculates zd%
      
      fld st0
      fmulp st1,st0         ;zd%*zd%
      fld qword [yd%]
      fld st0
      fmulp st1,st0         ;yd%*yd%
      fld qword [xd%]
      fld st0
      fmulp st1,st0         ;xd%*xd%
      faddp st1,st0
      faddp st1,st0         ;npl%*npl%
      fsqrt
      fstp qword [npl%]
      
      mov eax,[Ref8%]
      fld qword Fnow%[eax]
      fild dword [vt%]
      fsubp st1,st0
      fild dword [hdis%]
      fmul st1,st0
      fild dword [y%]
      faddp st1,st0
      fild dword [vt%]
      fmulp st1,st0
      faddp st1,st0
      fstp qword [num%]
      
      ret
      ]
      endproc
      
      def proc_dmpa
      Ref8%=8*(!x% + (!dm%+1)*!y%) : rem Done
      indRef8%=Fnow% + Ref8%
      
      |xd%=!vt%*(!y%+!hdis%)+!hdis%*(|indRef8%-!vt%)
      |yd%=-!vt%*(!x%-(!dm%/2))
      |zd%=-!hdis%*(!x%-(!dm%/2))
      
      |npl%=sqr(|xd%*|xd%+|yd%*|yd%+|zd%*|zd%)
      
      |num%=!vt%*(!y%+!hdis%)+!hdis%*(|indRef8%-!vt%)
      endproc
      
      
      rem Assembly Routine for Alpha Calc.
      def proc_beth(opt%)
      P%=beth%
      
      [opt opt%
      
      finit
      fild dword [dm%]
      fld1
      fld1
      faddp st1,st0
      fdivp st1,st0
      fild dword [x%]
      fsubp st1,st0         ;dm%/2 - x% calculated here
      fchs
      fstp qword [xd%]
      
      fild dword [hdis%]
      fild dword [y%]
      faddp st1,st0
      fstp qword [yd%]
      
      mov eax,[Ref8%]
      fld qword Fnow%[eax]
      fild dword [vt%]
      fsubp st1,st0
      fst qword [zd%]
      fld st0
      fmulp st1,st0         ;zd%*zd%
      fld qword [yd%]
      fld st0
      fmulp st1,st0         ;yd%*yd%
      fld qword [xd%]
      fld st0
      fmulp st1,st0         ;xd%*xd%
      faddp st1,st0
      faddp st1,st0
      fsqrt
      fstp qword [rdis%]    ;rdis%
      
      fild dword [y%]
      fild dword [hdis%]
      fadd st1,st0
      fmulp st1,st0
      fstp qword [yfact%]
      
      mov eax,[Ref8%]
      fld qword Fnow%[eax]
      fild dword [vt%]
      fsubr st1,st0
      fmulp st1,st0
      fstp qword [zfact%]
      
      ret
      ]
      endproc
      
      def proc_dmpb
      Ref8%=8*(!x% + (!dm%+1)*!y%) : rem Done
      indRef8%=Fnow% + Ref8%
      
      |xd%=(!x%-(!dm%/2)) : |yd%=!y%+!hdis% : |zd%=|indRef8%-!vt%
      |rdis%=sqr(|xd%*|xd%+|yd%*|yd%+|zd%*|zd%)
      
      |yfact%=!hdis%*(!hdis%+!y%) : |zfact%=!vt%*(!vt%-|indRef8%)
      endproc