rem Field ..... Rev 5.3 rem A J Tooth / 25th July 2003 Modified 10th August *FLOAT64 rem Call Utility Procedures proc_fullscreen(22) : rem Set up use of Full Screen himem=lomem + 30000000 rem Identifies the Program path, for later fetching a .jpg image proc_getpath rem Place Origin and Draw Axes input " Vertical Height (100)",vt if vt=0 then vt=100 input " Horizontal Distance (350)",hdis if hdis = 0 then hdis=300 input " High or Low Resolution? ( -h- or -l-)(l)",in$ if in$="" then in$="l" input " How many smoothing passes (7)",Ps% if Ps%=0 then Ps%=7 print " Press -Enter- to use default .JPG for background." print " Otherwise enter name of picture. " print : print " * MUST be in same file as this program. " print : print " * MUST be in format ' Picname.jpeg ' or ' Picname.jpg ', INCLUDING EXTENSION " print : input " Enter Picture Name .... ",Pname$ if Pname$="" then Pname$="Sky3.jpeg" Pname$= path$ + Pname$ rem Set up square field of side dm% (256 or 512) if in$="h" then dm%=512 else dm%=256 rem Set up Initial Field Arrays rem Dimension Field Array and Projection Array print tab(15,15);"STATUS: Dimensioning Arrays " dim Fnow(dm%,dm%,2) : dim alpbet(dm%,dm%,2) : dim pos(dm%,dm%,2) dim trac(25000,6) : dim Fsm(dm%,dm%) dim Norm(dm%,dm%,2,3) colour 7 rem Centres Graphics and sets scale xmax%=1024 : zmax%=768 scl%=5000 ofset%=0 : Zmax%=175 fct=.5*rnd(1) : rem Used in Perturbation Procedure origin xmax%,ofset% vht=vt : R=sqr(vht*vht+hdis*hdis) rem Randomize landscape mouse on 2 print tab(15,15);"STATUS: Randomising Landscape " proc_land(1,1) rem Smooth the landscape ps%=0 while ps%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.5*fn_cosPh(x%,vy%) when 1 : fac=0 : rgb%=fn_Refl(x%,vy%,vht,hdis) endcase else fac=0 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(r%*(1.3+fac)*(1-t%)) Gre%=int(g%*(1.3+fac)*(1-t%)) Blu%=int(b%*(1.3+fac)) endif vdu 19,10,16,Red%,Gre%,Blu% gcol 0,10 plot 85,scl%*pos(x%,vy%,1),scl%*pos(x%,vy%,2) vy%=vy%-1 plot 85,scl%*pos(x%,vy%,1),scl%*pos(x%,vy%,2) vy%=vy%+1 next x% next y% sys "GetCurrentProcess" to hprocess% sys "SetPriorityClass", hprocess%, &20 endproc rem -------------------------------------------------- rem FUNCTIONS rem ========= rem Calculates line of sight distance for 2-D View Angle def fn_rdis(Vht,Hdis,Ex%,Ey%) local xd,yd,zd xd=(Ex%-(dm%/2)) : yd=Ey%+Hdis : zd=Fnow(Ex%,Ey%,0)-Vht =sqr(xd*xd+yd*yd+zd*zd) rem Calculates normalising factor for normal to position vector plane def fn_npl(Vht,Hdis,Ex%,Ey%) local xd,yd,zd xd=Vht*(Ey%+Hdis)+Hdis*(Fnow(Ex%,Ey%,0)-Vht) yd=-Vht*(Ex%-(dm%/2)) zd=-Hdis*(Ex%-(dm%/2)) =sqr(xd*xd+yd*yd+zd*zd) rem ============================================================ rem Calculates Cosine of Light Incidence to get Intensity def fn_cosPh(Xx%,Yy%) local Term1,Term2,Term3,k%,hit%,Tes1%,Tes2%,Tes3,TesH rem Check if the sunlight hits an intervening hill k%=0 : hit%=0 repeat k%+=5 rem Calculate test vector point along line of sunlight Tes1%=int(Xx% - k%*Svec1) Tes2%=int(Yy% - k%*Svec2) Tes3=(Fnow(Xx%,Yy%,0) - k%*Svec3) if (Tes1%>-1 and Tes1%<(dm%+1) and Tes2%>-1 and Tes2%<(dm%+1)) then TesH=Fnow(Tes1%,Tes2%,0) if Tes3=175 or hit%=1) Term1=Svec1*Norm(Xx%-1,Yy%,1,1) Term2=Svec2*Norm(Xx%-1,Yy%,1,2) Term3=Svec3*Norm(Xx%-1,Yy%,1,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 then ResCos=0 =ResCos 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 2480 proc_alpha(vht,hdis,x%,y%) proc_beta(vht,hdis,x%,y%) 2480 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%,0)-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=fn_rdis(Vht,Hdis,Xx%,Yy%) yfact=Hdis*(Hdis+Yy%) : zfact=Vht*(Vht-Fnow(Xx%,Yy%,0)) alpbet(Xx%,Yy%,2)=acs((yfact+zfact)/(rdis*R)) rem Convert from Polar to x/y Coords for Printing to Screen pos(Xx%,Yy%,1)=alpbet(Xx%,Yy%,2)*-sin(alpbet(Xx%,Yy%,1)) pos(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(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 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 Gets the Program path, so that the linked .jpg file can be found def proc_getpath local len% dim path% 255 sys "GetModuleFileName", 0, path%, 255 to len% repeat len% -= 1 : until path%?len% = asc"\" path%?(len%+1) = 13 path$ = $path% endproc rem +++++++++++++++++++++++++++++++++++++++++++++++++++ rem Randomizes landscape def proc_land(Lps%,Flg%) local nd%,ind%,Max%,frstx%,drsty%,cntx%,cnty%,scndx%,scndy% sys "GetCurrentProcess" to hprocess% sys "SetPriorityClass", hprocess%, &100 Max%=150 stp%=0 nd%=0 if Flg%=1 then rem Set end-points xl%=0 : yl%=0 : xh%=dm% : yh%=dm% proc_set(xl%,yl%,xh%,yh%) 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 if Fnow(cntx%,frsty%,1)=1 then scndx%=cntx% else cntx%+=1 until (cntx%=(dm%+1) or scndx%>0) repeat if Fnow(frstx%,cnty%,1)=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(frstx%,frsty%,scndx%,scndy%,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 3780 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 3780 until (nd%=1) sys "GetCurrentProcess" to hprocess% sys "SetPriorityClass", hprocess%, &20 endproc rem +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def proc_water local x%,y% for x%=0 to dm% for y%=0 to dm% if Fnow(x%,y%,0)<0 then Fnow(x%,y%,0)=-1 next y% Fnow(x%,dm%,0)=-1 next x% endproc rem +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rem Perturbs mid-point def proc_midl(Xl%,Yl%,Xh%,Yh%,Fl%) local fact,red%,avg,xmid%,ymid% rem Reduce max perturbation with each loop. fact=1.2+fct red%=(fact)^(1.1*lp%) xmid%=Xl%+((Xh%-Xl%)/2) : ymid%=Yl%+((Yh%-Yl%)/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%=sw%*rnd(int(Max%/red%)) Fnow(xmid%,ymid%,1)=1 avg=(Fnow(Xl%,Yl%,0)+Fnow(Xl%,Yh%,0)+Fnow(Xh%,Xl%,0)+Fnow(Xh%,Yh%,0))/4 Fnow(xmid%,ymid%,0)=int(avg+pert%) xnew%=xmid% : ynew%=ymid% endproc rem ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ rem Initial Corner Points def proc_set(Xl%,Yl%,Xh%,Yh%) Fnow(Xl%,Yl%,1)=1 : Fnow(Xh%,Yh%,1)=1 Fnow(Xl%,Yh%,1)=1 : Fnow(Xh%,Yl%,1)=1 endproc 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 Ext%=6*Fnow(xn%,yn%,0) Fnow(xn%,yf%,1)=1 : Fnow(xn%,yf%,0)=(Ext%+Fnow(xf%,yf%,0)+Fnow(xs%,yf%,0))/8 Fnow(xf%,yn%,1)=1 : Fnow(xf%,yn%,0)=(Ext%+Fnow(xf%,yf%,0)+Fnow(xf%,ys%,0))/8 Fnow(xn%,ys%,1)=1 : Fnow(xn%,ys%,0)=(Ext%+Fnow(xf%,ys%,0)+Fnow(xs%,ys%,0))/8 Fnow(xs%,yn%,1)=1 : Fnow(xs%,yn%,0)=(Ext%+Fnow(xs%,yf%,0)+Fnow(xs%,ys%,0))/8 endproc rem -------------------------------------------------------------- rem Smooths the Field def proc_smooth local x%,y%,a%,xx%,yy% sys "GetCurrentProcess" to hprocess% sys "SetPriorityClass", hprocess%, &100 for a%=0 to dm% Fsm(0,a%)=0 : Fsm(dm%,a%)=0 : Fsm(a%,0)=0 : Fsm(a%,dm%)=0 next a% for x%=1 to (dm%-1) for y%=1 to (dm%-1) Fsm(x%,y%)=Fnow(x%-1,y%-1,0)+Fnow(x%-1,y%,0)+Fnow(x%-1,y%+1,0) Fsm(x%,y%)=Fsm(x%,y%)+Fnow(x%,y%-1,0)+Fnow(x%,y%,0)+Fnow(x%,y%+1,0) Fsm(x%,y%)=Fsm(x%,y%)+Fnow(x%+1,y%-1,0)+Fnow(x%+1,y%,0)+Fnow(x%+1,y%+1,0) Fsm(x%,y%)=Fsm(x%,y%)/9 next y% next x% for xx%=0 to dm% for yy%=0 to dm% Fnow(xx%,yy%,0)=Fsm(xx%,yy%) next yy% next xx% sys "GetCurrentProcess" to hprocess% sys "SetPriorityClass", hprocess%, &20 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 local x%,y%,m1,m2 sys "GetCurrentProcess" to hprocess% sys "SetPriorityClass", hprocess%, &100 for y%=dm% to 1 step -1 for x%=0 to (dm%-1) p11=x% : p12=y% : p13=Fnow(x%,y%,0) p21=x% : p22=y%-1 : p23=Fnow(x%,y%-1,0) p31=x%+1 : p32=y% : p33=Fnow(x%+1,y%,0) p41=x%+1 : p42=y%+1 : p43=Fnow(x%+1,y%-1,0) q11=p21-p11 : q12=p22-p12 : q13=p23-p13 q21=p31-p21 : q22=p32-p22 : q23=p33-p23 q31=p41-p31 : q32=p42-p32 : q33=p43-p33 m1=sqr((q12*q23-q13*q22)^2 + (q13*q21-q11*q23)^2 + (q11*q22-q12*q21)^2) Norm(x%,y%,1,1)=-(q12*q23-q13*q22)/m1 Norm(x%,y%,1,2)=-(q13*q21-q11*q23)/m1 Norm(x%,y%,1,3)=-(q11*q22-q12*q21)/m1 next x% next y% sys "GetCurrentProcess" to hprocess% sys "SetPriorityClass", hprocess%, &20 endproc rem +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def proc_pic(Edge%,Pic$) local pas% pas%=0 picture$ = Pic$ : rem. Image filename 5260 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 lm% picture%!(2*lm%) = 0 iid%!0 = &7BF80980 iid%!4 = &101ABF32 iid%!8 = &AA00BB8B iid%!12 = &AB0C3000 sys olpp%, picture%, 0, 0, 0, iid%, gpp% pas%+=1 if (pas%=1 and !gpp% = 0) then lg%=len(picture$) picture$=left$(picture$,(lg%-5)) + ".jpg" endif if (pas%=1 and !gpp%=0) goto 5260 if (pas%=2 and !gpp%=0) then error 0, "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) local Rgb%,k%,lam,hit% Dx%=Xx%-(dm%/2) rem (q1,q2,q3) is the line-of-sight vector to the point on the SEA ONLY. q1=Dx% : q2=Yy%+Hdis : q3=-Vht 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=(Yy%+Hdis)/Mq : s3=Vht/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%-Yy%)/(Yy%+Hdis) lam=(dm%-Yy%)*Mq/(Yy%+Hdis) rem Does the reflection hit intervening land, or carry on to meet the sky? k%=0 repeat k%+=5 Tes1%=Xx%+int(k%*s1) : Tes2%=Yy%+int(k%*s2) :Tes3=k%*s3 if (Tes1%>-1 and Tes1%<(dm%+1) and Tes2%>-1 and Tes2%<(dm%+1)) then if Fnow(Tes1%,Tes2%,0)>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(vht,hdis,Tes1%,Tes2%,com) proc_bet(vht,hdis,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(Vht,Hdis,Xx%,Yy%,Com) local moda,num,xd,yd,zd rem Calculates modulus of vector from origin to rem where the reflection "hits" the sky. xd=(-Hdis*Vht*Com) - (dm%*Vht) yd=Dx%*Vht*(1+Com) zd=Dx%*Hdis*(1+Com) moda=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=(Vht*Hdis*Com) + (dm%*Vht) if moda=0 then alp=pi/2 else alp=acs(num/moda) if (Xx%-(dm%/2))>0 then alp=-alp endproc rem----------------------------------------------------- rem 2-D Projection View-Angle def proc_bet(Vht,Hdis,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=Vht*(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=Vht*Vht*(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