rem Reflection ..... Rev 2.4
      rem A J Tooth / 24th August 2004
      rem Revised 9th November 2006
      
      on error if (err=17) then quit
      
      himem=lomem + 100000000
      *FLOAT 64
      rem First choose a background picture from a list
      mode 22
      colour 132,0,0,50 : colour 4,0,0,50 : colour 132 : cls
      
      rem Choose a picture
      proc_pichoose
      
      proc_fullscreen(xscreen%,yscreen%) : rem Set up use of Full Screen
      
      rem Place Origin and Draw Axes
      proc_setup
      
      rem Sets surface pattern of WATER.
      proc_water(ch$,in$)
      
      rem Project waterscape plane onto viewing screen
      print tab(15,15);"STATUS: Calculating Perspective       "
      
      proc_pers(vht,hdis)
      
      print tab(15,15);"STATUS: Calculating Normals       "
      rem Calculates the normals to each of the surface elements
      proc_norm
      
      cls : mouse off
      
      rem Draw the waterscape
      repeat
        cls : proc_Fdraw
        
        rem Press any key to exit
        sp$=get$
        rem Permits changing view position
        if sp$="h" then vht=int(1.5*vht) : proc_pers(vht,hdis)
        if sp$="l" then vht=int(0.5*vht) : proc_pers(vht,hdis)
        if sp$="f" then hdis=int(0.75*hdis) : proc_pers(vht,hdis)
        if sp$="b" then hdis=int(1.25*hdis) : proc_pers(vht,hdis)
      until (sp$<>"h" and sp$<>"l" and sp$<>"f" and sp$<>"b")
      if sp$<>" " then run
      
      quit
      end
      rem End of Programme ---------------------------------------
      
      rem Draw Points on Screen
      def proc_Fdraw
      local x%,y%,a$,xhorL%,xhorR%,diff%,g%
      
      edge%=int(scl%*pns(0,dm%,2))
      rem Prints the above JPG file to the screen.
      proc_pic(edge%)
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &100
      
      rem Displays Terrain
      cld%=0
      for y%=dm% to 1 step -1
        vy%=y%
        move scl%*pns(0,vy%,1),scl%*pns(0,vy%,2)
        move scl%*pns(0,vy%-1,1),scl%*pns(0,vy%-1,2)
        for x%=1 to dm%
          rem Calculate the reflection
          ch=Fnow(x%,vy%)
          !work%=fn_Refl(x%,vy%,vht,hdis,ch)
          
          rem Set r/g/b levels
          if !work%=-1 then !work%=0
          Red&=?(work%)
          Gre&=?(work%+1)
          Blu&=?(work%+2)
          
          rem Set screen colour, and plot two adjacent triangles in that colour
          colour 10,Red&,Gre&,Blu& : gcol 10
          plot 85,scl%*pns(x%,vy%,1),scl%*pns(x%,vy%,2)
          vy%=vy%-1
          plot 85,scl%*pns(x%,vy%,1),scl%*pns(x%,vy%,2)
          vy%=vy%+1
        next x%
      next y%
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &20
      endproc
      rem --------------------------------------------------
      
      rem FUNCTIONS
      rem Calculates normalising factor for normal to position vector plane
      def fn_npl(Vht,Hdis,Ex%,Ey%)
      local xd,yd,zd,E
      
      xd=Vht*(Ey%+Hdis)+Hdis*(Fnow(Ex%,Ey%)-Vht)
      E=1.0*(Ex%-(dm%/2))
      yd=-Vht*E
      zd=-Hdis*E
      =sqr(xd*xd+yd*yd+zd*zd)
      rem ----------------------------------------------------------
      
      rem PROCEDURES
      rem ==========
      rem Calculates perspective
      def proc_pers(vht,hdis)
      local x%,y%
      print tab(15,15);"STATUS: Calculating Perspective       "
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &100
      
      R=sqr(vht*vht+hdis*hdis)
      for x%=0 to dm%
        for y%=0 to dm%
          skp%=0
          on error local proc_error
          if skp%=1 then goto 1610
          
          proc_alpha(vht,hdis,x%,y%)
          proc_beta(vht,hdis,x%,y%)
          
 1610   next y%
      next x%
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &20
      endproc
      rem-----------------------------------------------------
      
      rem Angle between vertical plane and plane containing position vector
      def proc_alpha(Vht,Hdis,Xx%,Yy%)
      local npl,num
      
      npl=fn_npl(Vht,Hdis,Xx%,Yy%)
      num=Vht*(Yy%+Hdis)+Hdis*(Fnow(Xx%,Yy%)-Vht)
      if npl=0 then alpbet(Xx%,Yy%,1)=pi/2 else alpbet(Xx%,Yy%,1)=acs(num/npl)
      
      if (Xx%-(dm%/2))>0 then alpbet(Xx%,Yy%,1)=-alpbet(Xx%,Yy%,1)
      endproc
      rem-----------------------------------------------------
      
      rem 2-D Projection View-Angle
      def proc_beta(Vht,Hdis,Xx%,Yy%)
      local rdis
      
      rdis=sqr((Xx%-(dm%/2))^2 + (Yy%+Hdis)^2 + (Fnow(Xx%,Yy%)-Vht)^2)
      yfact=Hdis*(Hdis+Yy%) : zfact=Vht*(Vht-Fnow(Xx%,Yy%))
      alpbet(Xx%,Yy%,2)=acs((yfact+zfact)/(rdis*R))
      
      rem Convert from Polar to x/y Coords for Printing to Screen
      pns(Xx%,Yy%,1)=alpbet(Xx%,Yy%,2)*-sin(alpbet(Xx%,Yy%,1))
      pns(Xx%,Yy%,2)=alpbet(Xx%,Yy%,2)*cos(alpbet(Xx%,Yy%,1))
      endproc
      rem --------------------------------------------------
      
      rem Utility Procedures
      rem ==================
      
      rem Set up use of Full Screen
      def proc_fullscreen(return xscreen%,return yscreen%)
      
      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 Define the water surface
      def proc_water(Ch$,In$)
      local x%,y%,p,q
      
      case Ch$ of
        when "c" :
          rem Case of CIRCULAR waves.
          case In$ of
            when "l" : fac%=25
            when "h" : fac%=40
            when "u" : fac%=75
            otherwise error
          endcase
          rem Approximation to a wave-packet. Sinusoidal, not Bessel
          for x%=0 to dm%
            for y%=0 to dm%
              p=1.0*(x%-dm%/2) : q=1.0*(y%-dm%/2)
              r=sqr(p*p + q*q)
              Fnow(x%,y%)=(exp(-(((r-50)/fac%)^2)))*Wv*cos(r/Fr)
            next y%
          next x%
        when "p" :
          rem Case of PLANE waves
          for x%=0 to dm%
            for y%=0 to dm%
              Fnow(x%,y%)=Wv*cos(y%/Fr)
            next y%
          next x%
        when "n" :
          rem Case of NO waves
          for x%=0 to dm%
            for y%=0 to dm%
              Fnow(x%,y%)=0
            next y%
          next x%
      endcase
      
      endproc
      rem +++++++++++++++++++++++++++++++++++++++++++++++++++
      
      rem Calculate Normals to every triangulation
      def proc_norm
      local x%,y%,m,p1,p2,p3
      
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &100
      
      for y%=dm% to 1 step -1
        for x%=0 to (dm%-1)
          
          q13=Fnow(x%,y%-1)-Fnow(x%,y%)
          q23=Fnow(x%+1,y%)-Fnow(x%,y%-1)
          
          p1=-q23-q13
          
          m=sqr(p1*p1 + q13*q13 + 1.0)
          
          rem Normal UNIT vector
          Norm(x%,y%,1)=p1/m
          Norm(x%,y%,2)=q13/m
          Norm(x%,y%,3)=1.0/m
          
        next x%
      next y%
      sys "GetCurrentProcess" to hprocess%
      sys "SetPriorityClass", hprocess%, &20
      endproc
      rem +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      
      rem Loads a .jpg picture for the sky backdrop
      def proc_pic(Edge%)
      
      picture$ =pre$ + pic$        : rem Image filename
      xpos% = 0                    : rem X position (pixels)
      ypos% = 0                    : rem Y position (pixels)
      xsize% = 1024                : rem Width (pixels)
      ysize% = 768-((Edge%/2)-10)  : 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 lmn%
      picture%!(2*lmn%) = 0
      
      iid%!0  = &7BF80980
      iid%!4  = &101ABF32
      iid%!8  = &AA00BB8B
      iid%!12 = &AB0C3000
      
      sys olpp%, picture%, 0, 0, 0, iid%, gpp%
      if !gpp% = 0 error 0, "OleLoadPicturePath failed"
      
      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(Xx%,Yy%,Vht,Hdis,Ch)
      local Rgb%,lam,q1,q2,q3,resn,Tes1%,Tes2%,Tes3
      local s1,s2,s3,Ms,Dx%
      
      Dx%=Xx%-(dm%/2)
      
      rem (q1,q2,q3) is the line-of-sight vector to the point on the SEA ONLY.
      q1=1.0*Dx% : q2=1.0*(Yy%+Hdis) : q3=1.0*(Ch-Vht)
      resn=q1*Norm(Xx%,Yy%,1) + q2*Norm(Xx%,Yy%,2) + q3*Norm(Xx%,Yy%,3)
      
      rem (s1,s2,s3) is the reflected vector FROM the point on the sea.
      s1=q1-(2*resn*Norm(Xx%,Yy%,1))
      s2=q2-(2*resn*Norm(Xx%,Yy%,2))
      s3=q3-(2*resn*Norm(Xx%,Yy%,3))
      
      rem Normalises to get the UNIT vector
      Ms=sqr(s1*s1 + s2*s2 + s3*s3)
      s1=s1/Ms : s2=s2/Ms : s3=s3/Ms
      
      rem lam is the factor determining where the reflection vector meets the sky
      lam=(dm%-Yy%)*Ms/(Yy%+Hdis-(2*resn*Norm(Xx%,Yy%,2)))
      
      Tes1%=Dx%+int(lam*s1) : Tes2%=Yy%+int(lam*s2) : Tes3=Ch + lam*s3
      
      rem Work out relevant angles
      proc_alp(vht,hdis,Tes1%,Tes2%,Tes3)
      
      rem Find out what colour the sky is at the point where
      rem the reflection vector "hits" the sky on screen
      =tint(posx,posy)
      rem-----------------------------------------------------
      
      rem Modified Functional Procedures specific to Reflection Routine
      rem =============================================================
      rem Angle between vertical plane and plane containing position vector
      def proc_alp(Vht,Hdis,Xx%,Yy%,Zz)
      local mdpxa,num,xpa
      
      rem Calculates modulus of vector PxA
      xpa=-Hdis*Zz - Vht*Yy%
      mdpxa=sqr(xpa*xpa + ((Vht*Vht + Hdis*Hdis)*Xx%*Xx%))
      
      rem Works out the angle between the plane containing the line of sight
      rem vector and the vertical plane.
      if mdpxa=0 then alp=pi/2 else alp=acs((Hdis*Zz + Yy%*Vht)/mdpxa)
      
      if Xx%>0 then alp=-alp
      
      rem Get the azimuthal angle
      proc_bet(vht,hdis,Xx%,Yy%,Zz)
      
      endproc
      rem-----------------------------------------------------
      
      rem 2-D Projection View-Angle
      def proc_bet(Vht,Hdis,Xa,Ya,Za)
      local rdis,xd,yd,zd,yfact,zfact,Y,Z
      
      rem Calculates modulus of line-of-sight vector to
      rem where the reflection vector hits the sky.
      Y=Ya+Hdis : Z=Za-Vht
      rdis=sqr(Xa*Xa+Y*Y+Z*Z)
      
      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
      bet=acs((Hdis*Y-Vht*Z)/(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 Select .bmp file
      def proc_pichoose
      local g%,rn%,n$
      
      dim pf% 75, ff% 30, fm% 255
      !pf% = 76
      pf%!4 = @hwnd%
      pf%!12 = ff%
      pf%!28 = fm%
      pf%!32 = 256
      pf%!52 = 6
      
      $ff% ="Images"+chr$0+"*.bmp;*.gif;*.jpg"+chr$0+chr$0
      sys "GetOpenFileName", pf% to result%
      
      if result%<>0 then fullname$=$$fm%
      rn%=len(fullname$)
      g%=0 : pic$=""
      repeat
        n$=mid$(fullname$,rn%-g%,1)
        if n$<>"\" then pic$=n$+pic$ : g%+=1
      until n$="\"
      
      pre$=left$(fullname$,rn%-g%)
      rem Change the current Directory
      command$="CD "+chr$(34)+pre$+chr$(34)
      oscli command$
      endproc
      rem ==================================================================
      
      rem Place Origin and Draw Axes
      def proc_setup
      input " Vertical Height (450)",vt
      if vt=0 then vt=450
      input " Horizontal Distance (1300)",hdis
      if hdis=0 then hdis=1300
      input " Wave Height (0.25)",Wv
      if Wv=0 then Wv=0.25
      input " Wavelength (3)",Fr
      if Fr=0 then Fr=3
      input " Ultra-High, High or Low Resolution? ( -u- or -h- or -l-)(h)",in$
      if in$="" then in$="h"
      input " Plane Wave or Circular ( -p- or -c- or -n-)(n)",ch$
      if ch$="" then ch$="n"
      
      rem Set up square field of side dm% (256 or 512 or 1024)
      case in$ of
        when "l" : dm%=256
        when "h" : dm%=512
        when "u" : dm%=1024 : vt=int(vt/2)
        otherwise
          error
      endcase
      
      rem Set up Initial Field Arrays
      rem Dimension Field Array and Projection Array
      print tab(15,15);"STATUS: Dimensioning Arrays       "
      dim Fnow(dm%,dm%) : dim alpbet(dm%,dm%,2) : dim pns(dm%,dm%,2)
      dim Norm(dm%,dm%,3)
      dim work% 3
      
      colour 7
      rem Centres Graphics and sets scale : REM Sets various parameters
      xmax%=xscreen% : zmax%=yscreen%
      scl%=8000
      ffset%=0 : Zmax%=175
      vht=vt : R=sqr(vht*vht+hdis*hdis)
      
      origin xmax%,ffset%
      endproc
      rem ==================================================================