Planet generating program

This code has been declared by its author to be Public Domain code.

Planet generating program by markcw

2006

This is a translation of planet.c by Torben AE. Mogensen.
It generates planet maps in various projections.
It is a result of this thread http://www.blitzbasic.com/Community/posts.php?topic=63686
You may need the original to understand the options. http://www.diku.dk/~torbenm/

edit: updated this and replaced the commandline interface to just setting variables.
i also improved the search feature and made everything internal, so no color files, etc. needed.

;Translation of planet.c, planet generating program
;Copyright 1988-2003 Torben AE. Mogensen
;version of October 22 2003

;The program generates planet maps based on recursive spatial subdivision
;of a tetrahedron containing the globe. The output is BMP or PPM.

;The colours may optionally be modified according to latitude to move the
;icecaps lower closer to the poles, with a corresponding change in land
;colours.
;The Mercator map at magnification 1 is scaled to fit the width
;it uses the full height (it could extend infinitely)
;The orthographic projections are scaled so the full view would use the
;full Height. Areas outside the globe are coloured black.
;Stereographic and gnomic projections use the same scale as orthographic
;in the center of the picture, but distorts scale away from the center.

;It is assumed that pixels are square
;I have included procedures to print the maps as bmp or ppm
;bitmaps (portable pixel map) on standard output or specified files.

;I have tried to avoid using machine specific features, so it should
;be easy to port the program to any machine. Beware, though that due
;to different precision on different machines, the same seed numbers
;can yield very different planets.
;The primitive user interface is a result of portability concerns

Graphics 640,480,0,2
SetBuffer BackBuffer()

;these values can be changed to change world characteristics

Global debug = 1 ;show debug info
Global colorscheme = 1 ;specify colour scheme 0..5, default = 0
Global altColors = 0 ;Switch to alternative colour scheme

Global doshade = 0 ;Use 'bumpmap' shading
Global shade_angle# = 0 ;Angle of 'light' in bumpmap shading, eg. 150.0

Global Width = 512 ;Width in pixels, default = 256
Global Height = 256 ;Height in pixels, default = 128

Global view$ = "q" ;Projection: m=mercator, p=peter, q=squarep,
 ;M=Mollweide, S=Sinusoid, s=stereographic, o=orthographic
 ;g=gnomonic, a=azimuth, c=conical, h=heightfield, f=search

Global rseed# = 0.2 ;Seed as number between 0.0..1.0, eg. 0.3
Global M# = 0.19 ;Initial altitude, default = 0.3
Global dd1# = 0.15 ;Weight for altitude difference, default = 0.4
Global dd2# = 0.03 ;Weight for distance, default = 0.03
Global POW# = 0.2 ;Power for distance function, default = 0.47

Global scale# = 1.0 ;Magnification, default = 1.0
Global longi# = 0.0 ;Longitude of centre in degrees, default = 0.0
Global lat# = 0.0 ;Latitude of centre in degrees, default = 0.0
Global vgrid# = 0.0 ;Vertical gridsize in degrees, 0.0 = no grid
Global hgrid# = 0.0 ;Horisontal gridsize in degrees, 0.0 = no grid

Global latic = 0 ;Colour depends on latitude, default = only altitude
Global lighter = 0 ;Use lighter colours (original scheme only)

Global do_outline = 0 ;Trace the edges of land in black on colour map
Global do_bw = 0 ;Produce a black and white outline map

Global filename$ = "planet" ;map image file
Global file_type = 0 ;map file types, 0=BMP/1=PPM

Global searchfile$ = "search.txt" ;map search text file
Global increment# = 0.01 ;seed search increment, eg. 0.00000001

;globals used in calculations
Global BLACK=0
Global WHITE=1
Global BLUE0=2
Global BLUE1, LAND0, LAND1, LAND2, LAND4;
Global GREEN1, BROWN0, GREY0;
Global BACK = WHITE ;BLACK

Global shade;
Global Depth ;depth of subdivisions
Global r1#,r2#,r3#,r4# ;seeds
Global best = 500000

Global cla#,sla#,clo#,slo#;

Global DEG2RAD#=0.0174532918661 ;Pi/180
Global Rad2Deg#=(180/Float(Pi))

;arrays used in calculations
Dim col(256,128)
Dim shades(256,128)
Dim heights(256,128)
Dim rtable(256)
Dim gtable(256)
Dim btable(256)
Dim outx(256*128)
Dim outy(256*128)
Dim cl0(60,30);
Dim weight(30);
Dim colors(9,3)
Dim moll_table#(32)

;set colors array from data
Select colorscheme
 Case 1: Restore atlas_colors
 Case 2: Restore mars_colors
 Case 3: Restore grey_colors
 Case 4: Restore moon_colors
 Case 5: Restore cyan_colors
 Default: Restore colors
End Select
For i=0 To 8 ;copy the 9 base colors
 For j=0 To 2: Read colors(i,j): Next
Next

Restore moll_data
For i=0 To 18 ;copy the 19 moll values
 Read moll_table#(i)
Next

If (view = "h")
 Dim heights(Width,Height) ;heightmap
Else
 Dim col(Width,Height) ;map colors
EndIf
If (doshade)
 Dim shades(Width,Height) ;map shading
EndIf

If (debug)
 Print "> -s "+rseed+" -i "+M+" -a "+dd1+" -d "+dd2+" -p "+POW
 Print "> -w "+Width+" -h "+Height+" -B "+doshade
EndIf

;print error if unknown
Select (view)
 Case "m" ;mercator
 Case "p" ;peter
 Case "q" ;squarep
 Case "M" ;Mollweide
 Case "S" ;Sinusoid
 Case "s" ;stereographic
 Case "o" ;orthographic
 Case "g" ;gnomonic
 Case "a" ;azimuth
 Case "c" ;conical
 Case "h" ;heightfield
 Case "f" ;search
 Default : Print "Unknown projection: "+view : print_error()
End Select

If (view = "c")
 If (lat = 0) view = "m";
 ;Conical approaches mercator when lat -> 0
 If (Abs(lat) >= Pi - 0.000001) view = "s";
 ;Conical approaches stereo when lat -> +/- 90
EndIf

;init calculations
If (longi>180) longi = longi - 360;
longi = longi*DEG2RAD;
lat = lat*DEG2RAD;

sla = Sin(lat): cla = Cos(lat);
slo = Sin(longi): clo = Cos(longi);

Depth = 3*(Int(log_2(scale*Height)))+6;
r1 = rseed;
r1 = rand2(r1,r1);
r2 = rand2(r1,r1);
r3 = rand2(r1,r2);
r4 = rand2(r2,r3);

setcolours() ;set rgb arrays

If (debug) Print "+----+----+----+----+----+"

Select (view)
 Case "m" : mercator() ;Mercator projection
 Case "p" : peter() ;Peters projection (area preserving cylindrical)
 Case "q" : squarep() ;Square projection (equidistant latitudes)
 Case "M" : mollweide() ;Mollweide projection (area preserving)
 Case "S" : sinusoid() ;Sinusoid projection (area preserving)
 Case "s" : stereo() ;Stereographic projection
 Case "o" : orthographic() ;Orthographic projection
 Case "g" : gnomonic() ;Gnomonic projection
 Case "a" : azimuth() ;Area preserving azimuthal projection
 Case "c" : conical() ;Conical projection (conformal)
 Case "h" : heightfield() ;heightfield
 Case "f" ;Search
 file=WriteFile(searchfile)
 While (rseed<1)
  Restore map0 ;set map data
  readmap()
  search(file)
  best = 500000 ;reset
  rseed = rseed + increment ;next seed
  r1 = rseed;
  r1 = rand2(r1,r1);
  r2 = rand2(r1,r1);
  r3 = rand2(r1,r2);
  r4 = rand2(r2,r3);
 Wend
 CloseFile file
End Select

If (do_outline) makeoutline(do_bw) ;outlines

If (doshade) smoothshades() ;bumpmap shading

;add filetype extension
Select (file_type)
Case 1 : filename=filename+".ppm" ;ppm
Default : filename=filename+".bmp" ;bmp
End Select

;save map
Select (file_type)
 Case 1 ;ppm
  If (do_bw)
   printppmBW(filename);
  ElseIf (view <> "h")
   printppm(filename);
  Else
   printheights(filename);
  EndIf
 Default ;bmp
  If (do_bw)
   printbmpBW(filename);
  ElseIf (view <> "h")
   printbmp(filename);
  Else
   printheights(filename);
  EndIf
End Select

;draw map
image=CreateImage(Width,Height)
LockBuffer ImageBuffer(image)
For y=0 To ImageHeight(image)-1
 For x=0 To ImageWidth(image)-1
  If view="h"
   rgb=heights(x,y)
  Else
   r=rtable(col(x,y)) Shl 16
   g=gtable(col(x,y)) Shl 8
   b=btable(col(x,y))
   rgb=r+g+b
   If doshade
    s = (shades(x,y)/3)+85 ;smooth bumpmap
    r = s*rtable(col(x,y))/150;
    If (r>255) r=255;
    g = s*gtable(col(x,y))/150;
    If (g>255) g=255;
    b = s*btable(col(x,y))/150;
    If (b>255) b=255;
    rgb=(r Shl 16)+(g Shl 8)+b
   EndIf
  EndIf
  WritePixelFast x,y,rgb,ImageBuffer(image)
 Next
Next
UnlockBuffer ImageBuffer(image)
SetBuffer BackBuffer()

DrawImage image,50,200

WaitKey
End

;data labels

.colors ;colorscheme 0
Data 0,0,255 ;Dark blue depths
Data 0,128,255 ;Light blue shores
Data 0,255,0 ;Light green lowlands
Data 64,192,16 ;Dark green highlands
Data 64,192,16 ;Dark green Mountains
Data 128,128,32 ;Brown stoney peaks
Data 255,255,255 ;White - peaks
Data 0,0,0 ;Black - Space
Data 0,0,0 ;Black - Lines

.atlas_colors ;colorscheme 1
Data 0,0,192 ;Dark blue depths
Data 0,128,255 ;Light blue shores
Data 0,96,0 ;Dark green Lowlands
Data 0,224,0 ;Light green Highlands
Data 128,176,0 ;Brown mountainsides
Data 128,128,128 ;Grey stoney peaks
Data 255,255,255 ;White - peaks
Data 0,0,0 ;Black - Space
Data 0,0,0 ;Black - Lines

.mars_colors ;colorscheme 2
Data 55,28,75 ;0
Data 136,92,80 ;1
Data 148,100,82 ;2
Data 232,164,90 ;3
Data 248,180,92 ;4
Data 165,128,108 ;5
Data 255,255,255 ;6
Data 0,0,0 ;7
Data 0,0,0 ;8

.grey_colors ;colorscheme 3
Data 96,96,96 ;0
Data 152,152,152 ;2
Data 144,144,144 ;1
Data 168,168,168 ;3
Data 160,160,160 ;4
Data 144,144,144 ;5
Data 128,128,128 ;6
Data 0,0,0 ;7
Data 0,0,0 ;8

.moon_colors ;colorscheme 4
Data 112,112,98 ;0
Data 168,168,147 ;2
Data 160,160,140 ;1
Data 184,184,161 ;3
Data 176,176,154 ;4
Data 160,160,140 ;5
Data 144,144,126 ;6
Data 0,0,0 ;7
Data 0,0,0 ;8

.cyan_colors ;colorscheme 5
Data 74,112,150 ;0
Data 126,168,210 ;2
Data 120,160,200 ;1
Data 138,184,230 ;3
Data 132,176,220 ;4
Data 120,160,200 ;5
Data 108,144,180 ;6
Data 0,0,0 ;7
Data 0,0,0 ;8

.moll_data ;mollweide gridlines
Data 0.0, 0.0685055811, 0.1368109534
Data 0.2047150027, 0.2720147303, 0.3385041213
Data 0.4039727534, 0.4682040106, 0.5309726991
Data 0.5920417499, 0.6511575166, 0.7080428038
Data 0.7623860881, 0.8138239166, 0.8619100185
Data 0.9060553621, 0.9453925506, 0.9783738403, 1.0;

.map0 ;map data resembling earth
Data "::;--*;;@O;.;..-:**:;;::"
Data "-*o@O;*;-,;;o*O@@@@@O--;"
Data "..,O@Oo-,.;*--o*@@@O*;,."
Data "...;O-;...;O**oo@@@*;..."
Data "....;*;-.,*@@O*;*;o-...."
Data ".....:oO;,.;O@;.:,---;,."
Data ";;....*@O,..oO-...;;*;;;"
Data ",.....;O:...;*:...,Ooo,,"
Data "......-;...........,,;;,"
Data ".......;................"
Data "oooooooOoooooooooooooooo"

;functions

Function min(x,y)
 If x<y Then Return x Else Return y
End Function

Function max(x,y)
 If x<y Then Return y Else Return x
End Function

Function fmin#(x#,y#)
 If x<y Then Return x Else Return y
End Function

Function fmax#(x#,y#)
 If x<y Then Return y Else Return x
End Function

Function setcolours()

 Local i;
 Local crow;
 Local nocols = 256;

 If (debug)
  For crow = 0 To 8 ;MAXCOL-1
   Write "c"+crow+"="+colors(crow,0)+","+colors(crow,1)+","+colors(crow,2)+" ";
   If Not (crow+1) Mod 3 Print " "
  Next
 EndIf

 If (altColors)

  If (nocols < 8) nocols = 8;

  ;This Color table tries To follow the coloring conventions of
  ;several atlases.
  ;
  ;The First two colors are reserved For black And white
  ;1/4 of the colors are blue For the sea, dark being deep
  ;3/4 of the colors are land, divided as follows:
  ;nearly 1/2 of the colors are greens, with the low being dark
  ;1/8 of the colors shade from green through brown To grey
  ;1/8 of the colors are shades of grey For the highest altitudes
  ;
  ;The minimum Color table is:
  ;0 Black
  ;1 White
  ;2 Blue
  ;3 Dark Green
  ;4 Green
  ;5 Light Green
  ;6 Brown
  ;7 Grey
  ;And doesn't look very good. Somewhere between 24 And 32 colors
  ;is where this scheme starts looking good. 256, of course, is best.
    
  LAND0 = max(nocols / 4, BLUE0 + 1);
  BLUE1 = LAND0 - 1;
  GREY0 = nocols - (nocols / 8);
  GREEN1 = min(LAND0 + (nocols / 2), GREY0 - 2);
  BROWN0 = (GREEN1 + GREY0) / 2;
  LAND1 = nocols - 1;

  rtable(BLACK) = colors(7,0);
  gtable(BLACK) = colors(7,0);
  btable(BLACK) = colors(7,0);
  rtable(WHITE) = colors(6,0);
  gtable(WHITE) = colors(6,1);
  btable(WHITE) = colors(6,2);
  rtable(BLUE0) = colors(0,0);
  gtable(BLUE0) = colors(0,1);
  btable(BLUE0) = colors(0,2);

  For i=BLUE0+1 To BLUE1;
   rtable(i) = (colors(0,0)*(BLUE1-i)+colors(1,0)*(i-BLUE0))/(BLUE1-BLUE0);
   gtable(i) = (colors(0,1)*(BLUE1-i)+colors(1,1)*(i-BLUE0))/(BLUE1-BLUE0);
   btable(i) = (colors(0,2)*(BLUE1-i)+colors(1,2)*(i-BLUE0))/(BLUE1-BLUE0);
  Next
  For i=LAND0 To GREEN1-1
   rtable(i) = (colors(2,0)*(GREEN1-i)+colors(3,0)*(i-LAND0))/(GREEN1-LAND0);
   gtable(i) = (colors(2,1)*(GREEN1-i)+colors(3,1)*(i-LAND0))/(GREEN1-LAND0);
   btable(i) = (colors(2,2)*(GREEN1-i)+colors(3,2)*(i-LAND0))/(GREEN1-LAND0);
  Next
  For i=GREEN1 To BROWN0-1
   rtable(i) = (colors(3,0)*(BROWN0-i)+colors(4,0)*(i-GREEN1))/(BROWN0-GREEN1);
   gtable(i) = (colors(3,1)*(BROWN0-i)+colors(4,1)*(i-GREEN1))/(BROWN0-GREEN1);
   btable(i) = (colors(3,2)*(BROWN0-i)+colors(4,2)*(i-GREEN1))/(BROWN0-GREEN1);
  Next
  For i=BROWN0 To GREY0-1
   rtable(i) = (colors(4,0)*(GREY0-i)+colors(5,0)*(i-BROWN0))/(GREY0-BROWN0);
   gtable(i) = (colors(4,1)*(GREY0-i)+colors(5,1)*(i-BROWN0))/(GREY0-BROWN0);
   btable(i) = (colors(4,2)*(GREY0-i)+colors(5,2)*(i-BROWN0))/(GREY0-BROWN0);
  Next
  For i=GREY0 To nocols-1
   rtable(i) = (colors(5,0)*(nocols-i)+(colors(6,0)+1)*(i-GREY0))/(nocols-GREY0);
   gtable(i) = (colors(5,1)*(nocols-i)+(colors(6,1)+1)*(i-GREY0))/(nocols-GREY0);
   btable(i) = (colors(5,2)*(nocols-i)+(colors(6,2)+1)*(i-GREY0))/(nocols-GREY0);
  Next

 Else

  rtable(BLACK) = 0;
  gtable(BLACK) = 0;
  btable(BLACK) = 0;
  rtable(WHITE) = 255;
  gtable(WHITE) = 255;
  btable(WHITE) = 255;

  Local r, c;
  Local x#;

  While (lighter>0)
   For r = 0 To 7-1
    For c = 0 To 3-1
     x = Sqr(Float(colors(r,c))/256.0);
     colors(r,c) = Int(240.0*x+16);
    Next
   Next
   lighter=lighter-1
  Wend

  BLUE1 = (nocols-4)/2+BLUE0;

  If (BLUE1=BLUE0)
   rtable(BLUE0) = colors(0,0);
   gtable(BLUE0) = colors(0,1);
   btable(BLUE0) = colors(0,2);
  Else
   For i=BLUE0 To BLUE1;
    rtable(i) = (colors(0,0)*(BLUE1-i)+colors(1,0)*(i-BLUE0))/(BLUE1-BLUE0);
    gtable(i) = (colors(0,1)*(BLUE1-i)+colors(1,1)*(i-BLUE0))/(BLUE1-BLUE0);
    btable(i) = (colors(0,2)*(BLUE1-i)+colors(1,2)*(i-BLUE0))/(BLUE1-BLUE0);
   Next
  EndIf
  LAND0 = BLUE1+1: LAND2 = nocols-2: LAND1 = (LAND0+LAND2+1)/2;
  For i=LAND0 To LAND1-1
   rtable(i) = (colors(2,0)*(LAND1-i)+colors(3,0)*(i-LAND0))/(LAND1-LAND0);
   gtable(i) = (colors(2,1)*(LAND1-i)+colors(3,1)*(i-LAND0))/(LAND1-LAND0);
   btable(i) = (colors(2,2)*(LAND1-i)+colors(3,2)*(i-LAND0))/(LAND1-LAND0);
  Next

  If (LAND1=LAND2)
   rtable(LAND1) = colors(4,0);
   gtable(LAND1) = colors(4,1);
   btable(LAND1) = colors(4,2);
  Else
   For i=LAND1 To LAND2;
    rtable(i) = (colors(4,0)*(LAND2-i)+colors(5,0)*(i-LAND1))/(LAND2-LAND1);
    gtable(i) = (colors(4,1)*(LAND2-i)+colors(5,1)*(i-LAND1))/(LAND2-LAND1);
    btable(i) = (colors(4,2)*(LAND2-i)+colors(5,2)*(i-LAND1))/(LAND2-LAND1);
   Next
  EndIf
  LAND4 = nocols-1;
  rtable(LAND4) = colors(6,0);
  gtable(LAND4) = colors(6,1);
  btable(LAND4) = colors(6,2);

 EndIf

End Function

Function makeoutline(do_bw)

 Local i,j,k;

 Dim outx(Width*Height)
 Dim outy(Width*Height)

 k=0;
 For i=1 To Width-2
  For j=1 To Height-2
   If ((col(i,j) >= BLUE0 And col(i,j) <= BLUE1) And (col(i-1,j) >= LAND0 Or col(i+1,j) >= LAND0 Or col(i,j-1) >= LAND0 Or col(i,j+1) >= LAND0 Or col(i-1,j-1) >= LAND0 Or col(i-1,j+1) >= LAND0 Or col(i+1,j-1) >= LAND0 Or col(i+1,j+1) >= LAND0))
    outx(k) = i: outy(k) = j: k=k+1 ;outy(k++)
   EndIf
  Next
 Next

 If (do_bw)
  For i=0 To Width-1
   For j=0 To Height-1
    If (col(i,j) <> BLACK) col(i,j) = WHITE;
    While (k>0) ;(k-- >0)
     col(outx(k),outy(k)) = BLACK
     k=k-1
    Wend
   Next
  Next
 EndIf

End Function

Function smoothshades()

 Local i,j;

 For i=0 To Width-3 ;i<Width-2
  For j=0 To Height-3;
   shades(i,j) = (4*shades(i,j)+2*shades(i,j+1)+2*shades(i+1,j)+shades(i+1,j+2)+4)/9;
  Next
 Next

End Function

Function mercator() ;-pm

 Local y#,scale1#,cos2#,theta1#,theta0#;
 Local i,j,k;

 y = Sin(lat*Rad2Deg) ;Sin(lat)
 y = (1.0+y)/(1.0-y);
 y = 0.5*Log(y);
 k = Int(0.5*y*Width*scale/Pi);

 For j = 0 To Height-1
  If (debug And ((j Mod (Height/25)) = 0)) Write view ;%=mod
  y = Pi*(2.0*(j-k)-Height)/Width/scale;
  y = Exp(2.0*y);
  y = (y-1.0)/(y+1.0);
  scale1 = scale*Width/Height/Sqr(1.0-y*y)/Pi;
  cos2 = Sqr(1.0-y*y);
  Depth = 3*(Int(log_2(scale1*Height)))+3;
  For i = 0 To Width-1
   theta1 = longi-0.5*Pi+Pi*(2.0*i-Width)/Width/scale;
   col(i,j) = planet0(Cos(theta1*Rad2Deg)*cos2,y,-Sin(theta1*Rad2Deg)*cos2);
   ;col(i,j) = planet0(Cos(theta1)*cos2,y,-Sin(theta1)*cos2);
   If (doshade) shades(i,j) = shade;
  Next
 Next

 If (hgrid <> 0.0) ;draw horisontal gridlines
  theta0 = 0.0
  While theta0>-90.0
   theta0=theta0-hgrid
   theta1 = theta0
   While theta1<90.0
    theta1=theta1+hgrid
    y = Sin(DEG2RAD*theta1*Rad2Deg) ;Sin(DEG2RAD*theta1)
    y = (1.0+y)/(1.0-y);
    y = 0.5*Log(y);
    j = Height/2+Int(0.5*y*Width*scale/Pi)+k;
    If (j>=0 And j<Height)
     For i = 0 To Width-1: col(i,j) = BLACK: Next;
    EndIf
   Wend
  Wend
 EndIf

 If (vgrid <> 0.0) ;draw vertical gridlines
  theta0 = 0.0
  While theta0>-360.0
   theta0=theta0-vgrid
   theta1 = theta0
   While theta1<360.0
    theta1=theta1+vgrid
    i = Int(0.5*Width*(1.0+scale*(DEG2RAD*theta1-longi)/Pi));
    If (i>=0 And i<Width)
     For j = 0 To Height-1: col(i,j) = BLACK: Next;
    EndIf
   Wend
  Wend
 EndIf

End Function

Function peter() ;-pp

 Local y#,cos2#,theta1#,scale1#,theta0#
 Local k,i,j,water,land

 y = 2.0*Sin(lat*Rad2Deg) ;2.0*Sin(lat)
 k = Int(0.5*y*Width*scale/Pi);
 water = 0: land = 0;

 For j = 0 To Height-1
  If (debug And ((j Mod (Height/25)) = 0)) Write view
  y = 0.5*Pi*(2.0*(j-k)-Height)/Width/scale;
  If (Abs(y)>1.0)
   For i = 0 To Width-1
    col(i,j) = BACK;
    If (doshade) shades(i,j) = 255;
   Next
  Else
   cos2 = Sqr(1.0-y*y);
   If (cos2>0.0)
    scale1 = scale*Width/Height/cos2/Pi;
    Depth = 3*(Int(log_2(scale1*Height)))+3;
    For i = 0 To Width-1
     theta1 = longi-0.5*Pi+Pi*(2.0*i-Width)/Width/scale;
     col(i,j) = planet0(Cos(theta1*Rad2Deg)*cos2,y,-Sin(theta1*Rad2Deg)*cos2);
     ;col(i,j) = planet0(Cos(theta1)*cos2,y,-Sin(theta1)*cos2);
     If (doshade) shades(i,j) = shade;
     If (col(i,j) < LAND0) Then water=water+1 Else land=land+1;
    Next
   EndIf
  EndIf
 Next

 If (debug)
  Print " ": Print "water percentage: "+(100*water/(water+land));
 EndIf

 If (hgrid <> 0.0) ;draw horisontal gridlines
  theta0 = 0.0
  While theta0>-90.0
  theta0=theta0-hgrid
  theta1 = theta0
   While theta1<90.0
    theta1=theta1+hgrid
    y = 2.0*Sin(DEG2RAD*theta1*Rad2Deg) ;2.0*Sin(DEG2RAD*theta1)
    j = Height/2+Int(0.5*y*Width*scale/Pi)+k;
    If (j>=0 And j<Height)
     For i = 0 To Width-1: col(i,j) = BLACK: Next;
    EndIf
   Wend
  Wend
 EndIf

 If (vgrid <> 0.0) ;draw vertical gridlines
  theta0 = 0.0
  While theta0>-360.0
   theta0=theta0-vgrid
   theta1 = theta0
   While theta1<360.0
    theta1=theta1+vgrid
    i = Int(0.5*Width*(1.0+scale*(DEG2RAD*theta1-longi)/Pi));
    If (i>=0 And i<Width)
     For j = max(0,Height/2-Int(Width*scale/Pi)+k) To min(Height,Height/2+Int(Width*scale/Pi)+k)-1
      col(i,j) = BLACK;
     Next
    EndIf
   Wend
  Wend
 EndIf

End Function

Function squarep() ;-pq

 Local y#,scale1#,theta1#,cos2#,theta0#;
 Local k,i,j;

 k = Int(lat*Width*scale/Pi);

 For j = 0 To Height-1
  If (debug And ((j Mod (Height/25)) = 0)) Write view
  y = (2.0*(j-k)-Height)/Width/scale*Pi;
  If (Abs(y)>0.5*Pi) ;(Abs(y)>=0.5*Pi)
   For i = 0 To Width-1
    col(i,j) = BACK;
    If (doshade) shades(i,j) = 255;
   Next
  Else
   cos2 = Cos(y*Rad2Deg) ;Cos(y)
   If (cos2>0.0)
    scale1 = scale*Width/Height/cos2/Pi;
    Depth = 3*(Int(log_2(scale1*Height)))+3;
    For i = 0 To Width-1
     theta1 = longi-0.5*Pi+Pi*(2.0*i-Width)/Width/scale;
     col(i,j) = planet0(Cos(theta1*Rad2Deg)*cos2,Sin(y*Rad2Deg),-Sin(theta1*Rad2Deg)*cos2);
     ;col(i,j) = planet0(Cos(theta1)*cos2,Sin(y),-Sin(theta1)*cos2);
     If (doshade) shades(i,j) = shade;
    Next
   EndIf
  EndIf
 Next

 If (hgrid <> 0.0) ;draw horisontal gridlines
  theta0 = 0.0
  While theta0>-90.0
  theta0=theta0-hgrid
  theta1 = theta0
   While theta1<90.0
    theta1=theta1+hgrid
    y = DEG2RAD*theta1;
    j = Height/2+Int(0.5*y*Width*scale/Pi)+k;
    If (j>=0 And j<Height)
     For i = 0 To Width-1: col(i,j) = BLACK: Next;
    EndIf
   Wend
  Wend
 EndIf

 If (vgrid <> 0.0) ;draw vertical gridlines
  theta0 = 0.0
  While theta0>-360.0
   theta0=theta0-vgrid
   theta1 = theta0
   While theta1<360.0
    theta1=theta1+vgrid
    i = Int(0.5*Width*(1.0+scale*(DEG2RAD*theta1-longi)/Pi));
    If (i>=0 And i<Width)
     For j = max(0,Height/2-Int(0.25*Pi*Width*scale/Pi)+k) To min(Height,Height/2+Int(0.25*Pi*Width*scale/Pi)+k)-1
      col(i,j) = BLACK;
     Next
    EndIf
   Wend
  Wend
 EndIf

End Function

Function mollweide() ;-pM

 Local x#,y#,y1#,zz#,scale1#,cos2#,theta1#,theta2#,theta0#;
 Local i,j,i1=1,k;

 For j = 0 To Height-1
  If (debug And ((j Mod (Height/25)) = 0)) Write view
  y1 = 2*(2.0*j-Height)/Width/scale;
  If (Abs(y1)>=1.0)
   For i = 0 To Width-1
    col(i,j) = BACK;
    If (doshade) shades(i,j) = 255;
   Next
  Else
   zz = Sqr(1.0-y1*y1);
   y = 2.0/Pi*(y1*zz+ASin(y1))
   cos2 = Sqr(1.0-y*y);
   If (cos2>0.0)
    scale1 = scale*Width/Height/cos2/Pi;
    Depth = 3*(Int(log_2(scale1*Height)))+3;
    For i = 0 To Width-1
     theta1 = Pi/zz*(2.0*i-Width)/Width/scale;
     If (Abs(theta1)>Pi)
      col(i,j) = BACK;
      If (doshade) shades(i,j) = 255;
     Else
      theta1 = theta1 + longi-0.5*Pi;
      col(i,j) = planet0(Cos(theta1*Rad2Deg)*cos2,y,-Sin(theta1*Rad2Deg)*cos2);
      ;col(i,j) = planet0(Cos(theta1)*cos2,y,-Sin(theta1)*cos2);
      If (doshade) shades(i,j) = shade;
     EndIf
    Next
   EndIf
  EndIf
 Next

 If (hgrid <> 0.0) ;draw horisontal gridlines
  theta0 = 0.0
  While theta0>-90.0
  theta0=theta0-hgrid
  theta1 = theta0
   While theta1<90.0
    theta1=theta1+hgrid
    theta2 = Abs(theta1);
    x = Floor(theta2/5.0): y = theta2/5.0-x;
    y = (1.0-y)*moll_table(Int(x))+y*moll_table(Int(x+1));
    If (theta1<0.0) y = -y;
     j = Height/2+Int(0.25*y*Width*scale);
     If (j>=0 And j<Height)
      For i = max(0,Width/2-Int(0.5*Width*scale*Sqr(1.0-y*y))) To min(Width,Width/2+Int(0.5*Width*scale*Sqr(1.0-y*y)))-1
      col(i,j) = BLACK;
      Next
     EndIf
   Wend
  Wend
 EndIf

 If (vgrid <> 0.0) ;draw vertical gridlines
  theta0 = 0.0
  While theta0>-360.0
   theta0=theta0-vgrid
   theta1 = theta0
   While theta1<360.0
    theta1=theta1+vgrid
     If (DEG2RAD*theta1-longi+0.5*Pi>-Pi And DEG2RAD*theta1-longi+0.5*Pi<=Pi)
      x = 0.5*(DEG2RAD*theta1-longi+0.5*Pi)*Width*scale/Pi;
      j = max(0,Height/2-Int(0.25*Width*scale));
      y = 2*(2.0*j-Height)/Width/scale;
      i = Width/2+x*Sqr(1.0-y*y);
      For j = max(0,Height/2-Int(0.25*Width*scale)) To min(Height,Height/2+Int(0.25*Width*scale))
       y1 = 2*(2.0*j-Height)/Width/scale;
       If (Abs(y1)<=1.0)
        i1 = Width/2+x*Sqr(1.0-y1*y1);
        If (i1>=0 And i1<Width) col(i1,j) = BLACK;
       EndIf
       If (Abs(y)<=1.0)
        If (i<i1)
         For k=i+1 To i1-1 ;for (k=i+1; k<i1; k++)
          If (k>0 And k<Width) col(k,j) = BLACK;
         Next
        ElseIf (i>i1)
         k=i-1
         While k>i1 ;for (k=i-1; k>i1; k--)
          If (k>=0 And k<Width) col(k,j) = BLACK;
          k=k-1
         Wend
        EndIf
       EndIf
       y = y1;
       i = i1;
      Next
     EndIf
    Wend
   Wend
  EndIf

End Function

Function sinusoid() ;-pS

 Local y#,theta1#,theta2#,cos2#,l1#,i1#,scale1#,theta0#;
 Local k,i,j,l,c;

 k = Int(lat*Width*scale/Pi);

 For j = 0 To Height-1
  If (debug And ((j Mod (Height/25)) = 0)) Write view
  y = (2.0*(j-k)-Height)/Width/scale*Pi;
  If (Abs(y)>=0.5*Pi)
   For i = 0 To Width-1
    col(i,j) = BACK;
    If (doshade) shades(i,j) = 255;
   Next
  Else
   cos2 = Cos(y*Rad2Deg) ;Cos(y)
   If (cos2>0.0)
    scale1 = scale*Width/Height/cos2/Pi;
    Depth = 3*(Int(log_2(scale1*Height)))+3;
    For i = 0 To Width-1
     l = i*12/Width;
     l1 = l*Width/12.0;
     i1 = i-l1;
     theta2 = longi-0.5*Pi+Pi*(2.0*l1-Width)/Width/scale;
     theta1 = (Pi*(2.0*i1-Width/12)/Width/scale)/cos2;
     If (Abs(theta1)>Pi/12.0)
      col(i,j) = BACK;
      If (doshade) shades(i,j) = 255;
     Else
      col(i,j) = planet0(Cos((theta1+theta2)*Rad2Deg)*cos2,Sin(y*Rad2Deg),-Sin((theta1+theta2)*Rad2Deg)*cos2);
      ;col(i,j) = planet0(Cos(theta1+theta2)*cos2,Sin(y),-Sin(theta1+theta2)*cos2);
      If (doshade) shades(i,j) = shade;
     EndIf
    Next
   EndIf
  EndIf
 Next

 If (hgrid <> 0.0) ;draw horisontal gridlines
  theta0 = 0.0
  While theta0>-90.0
   theta0=theta0-hgrid
   theta1 = theta0
   While theta1<90.0
    theta1=theta1+hgrid
    y = DEG2RAD*theta1;
    cos2 = Cos(y*Rad2Deg) ;Cos(y)
    j = Height/2+Int(0.5*y*Width*scale/Pi)+k;
    If (j>=0 And j<Height)
     For i = 0 To Width-1
      l = i*12/Width;
      l1 = l*Width/12.0;
      i1 = i-l1;
      theta2 = (Pi*(2.0*i1-Width/12)/Width/scale)/cos2;
      If (Abs(theta2)<=Pi/12.0) col(i,j) = BLACK;
     Next
    EndIf
   Wend
  Wend
 EndIf

 If (vgrid <> 0.0) ;draw vertical gridlines
  theta0 = 0.0
  While theta0>-360.0
   theta0=theta0-vgrid
   theta1 = theta0
   While theta1<360.0
    theta1=theta1+vgrid
    i = Int(0.5*Width*(1.0+scale*(DEG2RAD*theta1-longi)/Pi));
    If (i>=0 And i<Width)
     For j = max(0,Height/2-Int(0.25*Pi*Width*scale/Pi)+k) To min(Height,Height/2+Int(0.25*Pi*Width*scale/Pi)+k)-1
      y = (2.0*(j-k)-Height)/Width/scale*Pi;
      cos2 = Cos(y*Rad2Deg) ;Cos(y)
      l = i*12/Width;
      l1 = l*Width/12.0+Width/24.0;
      i1 = i-l1;
      c = l1+i1*cos2;
      If (c>=0 And c<Width) col(c,j) = BLACK;
     Next
    EndIf
   Wend
  Wend
 EndIf

End Function

Function stereo() ;-ps

 Local x#,y#,ymin#,ymax#,z#,zz#,x1#,y1#,z1#,theta1#,theta2#,theta0#;
 Local i,j;

 ymin = 2.0;
 ymax = -2.0;

 For j = 0 To Height-1
  If (debug And ((j Mod (Height/25)) = 0)) Write view
  For i = 0 To Width-1
   x = (2.0*i-Width)/Height/scale;
   y = (2.0*j-Height)/Height/scale;
   z = x*x+y*y;
   zz = 0.25*(4.0+z);
   x = x/zz;
   y = y/zz;
   z = (1.0-0.25*z)/zz;
   x1 = clo*x+slo*sla*y+slo*cla*z;
   y1 = cla*y-sla*z;
   z1 = -slo*x+clo*sla*y+clo*cla*z;
   If (y1 < ymin) ymin = y1;
   If (y1 > ymax) ymax = y1;
   col(i,j) = planet0(x1,y1,z1);
   If (doshade) shades(i,j) = shade;
  Next
 Next

 If (hgrid <> 0.0) ;draw horisontal gridlines
  theta0 = 0.0
  While theta0>-90.0
   theta0=theta0-hgrid
   theta1 = theta0
   While theta1<90.0
    theta1=theta1+hgrid
    y = Sin(DEG2RAD*theta1*Rad2Deg) ;Sin(DEG2RAD*theta1)
    If (ymin <= y And y <= ymax)
     zz = Sqr(1-y*y);
     theta2=-Pi
     While theta2<Pi
      theta2=theta2+0.5/Width/scale
      x = Sin(theta2*Rad2Deg)*zz ;Sin(theta2)*zz
      z = Cos(theta2*Rad2Deg)*zz ;Cos(theta2)*zz
      x1 = clo*x+slo*z;
      y1 = slo*sla*x+cla*y-clo*sla*z;
      z1 = -slo*cla*x+sla*y+clo*cla*z;
      If (Abs(z1)<1.0)
       i = 0.5*(Height*scale*2.0*x1*(1+z1)/(1.0-z1*z1)+Width);
       j = 0.5*(Height*scale*2.0*y1*(1+z1)/(1.0-z1*z1)+Height);
       If (i>=0 And i<Width And j>=0 And j<Height) col(i,j) = BLACK;
      EndIf
     Wend
    EndIf
   Wend
  Wend
 EndIf

 If (vgrid <> 0.0) ;draw vertical gridlines
  theta2=-Pi
  While theta2<Pi
   theta2=theta2+0.5/Width/scale
   y = Sin(theta2*Rad2Deg) ;Sin(theta2)
   If (ymin <= y And y <= ymax)
    theta1 = 0.0
    While theta1<360.0
     theta1=theta1+vgrid
     x = Sin(DEG2RAD*theta1*Rad2Deg)*Cos(theta2*Rad2Deg);
     z = Cos(DEG2RAD*theta1*Rad2Deg)*Cos(theta2*Rad2Deg);
     ;x = Sin(DEG2RAD*theta1)*Cos(theta2);
     ;z = Cos(DEG2RAD*theta1)*Cos(theta2);
     x1 = clo*x+slo*z;
     y1 = slo*sla*x+cla*y-clo*sla*z;
     z1 = -slo*cla*x+sla*y+clo*cla*z;
     If (Abs(z1)<1.0)
      i = 0.5*(Height*scale*2.0*x1*(1+z1)/(1-z1*z1)+Width);
      j = 0.5*(Height*scale*2.0*y1*(1+z1)/(1-z1*z1)+Height);
      If (i>=0 And i<Width And j>=0 And j<Height) col(i,j) = BLACK;
     EndIf
    Wend
   EndIf
  Wend
 EndIf

End Function

Function orthographic() ;-po

 Local x#,y#,z#,x1#,y1#,z1#,ymin#,ymax#,theta1#,theta2#,zz#,theta0#;
 Local i,j;

 ymin = 2.0;
 ymax = -2.0;

 For j = 0 To Height-1
  If (debug And ((j Mod (Height/25)) = 0)) Write view
  For i = 0 To Width-1
   x = (2.0*i-Width)/Height/scale;
   y = (2.0*j-Height)/Height/scale;
   If (x*x+y*y>1.0)
    col(i,j) = BACK;
    If (doshade) shades(i,j) = 255;
   Else
    z = Sqr(1.0-x*x-y*y);
    x1 = clo*x+slo*sla*y+slo*cla*z;
    y1 = cla*y-sla*z;
    z1 = -slo*x+clo*sla*y+clo*cla*z;
    If (y1 < ymin) ymin = y1;
    If (y1 > ymax) ymax = y1;
    col(i,j) = planet0(x1,y1,z1);
    If (doshade) shades(i,j) = shade;
   EndIf
  Next
 Next

 If (hgrid <> 0.0) ;draw horisontal gridlines
  theta0 = 0.0
  While theta0>-90.0
   theta0=theta0-hgrid
   theta1 = theta0
   While theta1<90.0
    theta1=theta1+hgrid
    y = Sin(DEG2RAD*theta1*Rad2Deg) ;Sin(DEG2RAD*theta1)
    If (ymin <= y And y <= ymax)
     zz = Sqr(1-y*y);
     theta2=-Pi
     While theta2<Pi
      theta2=theta2+0.5/Width/scale
      x = Sin(theta2*Rad2Deg)*zz ;Sin(theta2)*zz
      z = Cos(theta2*Rad2Deg)*zz ;Cos(theta2)*zz
      x1 = clo*x+slo*z;
      y1 = slo*sla*x+cla*y-clo*sla*z;
      z1 = -slo*cla*x+sla*y+clo*cla*z;
      If (z1<=0.0)
       i = 0.5*(Height*scale*x1+Width);
       j = 0.5*(Height*scale*y1+Height);
       If (i>=0 And i<Width And j>=0 And j<Height) col(i,j) = BLACK;
      EndIf
     Wend
    EndIf
   Wend
  Wend
 EndIf

 If (vgrid <> 0.0) ;draw vertical gridlines
  theta2=-Pi
  While theta2<Pi
   theta2=theta2+0.5/Width/scale
   y = Sin(theta2*Rad2Deg) ;Sin(theta2)
   If (ymin <= y And y <= ymax)
    theta1 = 0.0
    While theta1<360.0
     theta1=theta1+vgrid
     x = Sin(DEG2RAD*theta1*Rad2Deg)*Cos(theta2*Rad2Deg);
     z = Cos(DEG2RAD*theta1*Rad2Deg)*Cos(theta2*Rad2Deg);
     ;x = Sin(DEG2RAD*theta1)*Cos(theta2);
     ;z = Cos(DEG2RAD*theta1)*Cos(theta2);
     x1 = clo*x+slo*z;
     y1 = slo*sla*x+cla*y-clo*sla*z;
     z1 = -slo*cla*x+sla*y+clo*cla*z;
     If (0.0>=z1)
      i = 0.5*(Height*scale*x1+Width);
      j = 0.5*(Height*scale*y1+Height);
      If (i>=0 And i<Width And j>=0 And j<Height) col(i,j) = BLACK;
     EndIf
    Wend
   EndIf
  Wend
 EndIf

End Function

Function gnomonic() ;-pg

 Local x#,y#,z#,x1#,y1#,z1#,zz#,theta1#,theta2#,ymin#,ymax#,theta0#;
 Local i,j;

 ymin = 2.0;
 ymax = -2.0;

 For j = 0 To Height-1
  If (debug And ((j Mod (Height/25)) = 0)) Write view
  For i = 0 To Width-1
   x = (2.0*i-Width)/Height/scale;
   y = (2.0*j-Height)/Height/scale;
   zz = Sqr(1.0/(1.0+x*x+y*y));
   x = x*zz;
   y = y*zz;
   z = Sqr(1.0-x*x-y*y);
   x1 = clo*x+slo*sla*y+slo*cla*z;
   y1 = cla*y-sla*z;
   z1 = -slo*x+clo*sla*y+clo*cla*z;
   If (y1 < ymin) ymin = y1;
   If (y1 > ymax) ymax = y1;
   col(i,j) = planet0(x1,y1,z1);
   If (doshade) shades(i,j) = shade;
  Next
 Next

 If (hgrid <> 0.0) ;draw horisontal gridlines
  theta0 = 0.0
  While theta0>-90.0
   theta0=theta0-hgrid
   theta1 = theta0
   While theta1<90.0
    theta1=theta1+hgrid
    y = Sin(DEG2RAD*theta1*Rad2Deg) ;Sin(DEG2RAD*theta1)
    If (ymin <= y And y <= ymax)
     zz = Sqr(1-y*y);
     theta2=-Pi
     While theta2<Pi
      theta2=theta2+0.5/Width/scale
      x = Sin(theta2*Rad2Deg)*zz ;Sin(theta2)*zz
      z = Cos(theta2*Rad2Deg)*zz ;Cos(theta2)*zz
      x1 = clo*x-slo*z;
      y1 = slo*sla*x+cla*y+clo*sla*z;
      z1 = slo*cla*x-sla*y+clo*cla*z;
      If (z1<>0.0)
       i = 0.5*(Height*scale*x1/z1+Width);
       j = 0.5*(Height*scale*y1/z1+Height);
       If (i>=0 And i<Width And j>=0 And j<Height) col(i,j) = BLACK;
      EndIf
     Wend
    EndIf
   Wend
  Wend
 EndIf

 If (vgrid <> 0.0) ;draw vertical gridlines
  theta2=-Pi
  While theta2<Pi
   theta2=theta2+0.5/Width/scale
   y = Sin(theta2*Rad2Deg) ;Sin(theta2)
   If (ymin <= y And y <= ymax)
    theta1 = 0.0
    While theta1<360.0
     theta1=theta1+vgrid
     x = Sin(DEG2RAD*theta1*Rad2Deg)*Cos(theta2*Rad2Deg);
     z = Cos(DEG2RAD*theta1*Rad2Deg)*Cos(theta2*Rad2Deg);
     ;x = Sin(DEG2RAD*theta1)*Cos(theta2);
     ;z = Cos(DEG2RAD*theta1)*Cos(theta2);
     x1 = clo*x-slo*z;
     y1 = slo*sla*x+cla*y+clo*sla*z;
     z1 = slo*cla*x-sla*y+clo*cla*z;
     If (z1<>0.0)
      i = 0.5*(Height*scale*x1/z1+Width);
      j = 0.5*(Height*scale*y1/z1+Height);
      If (i>=0 And i<Width And j>=0 And j<Height) col(i,j) = BLACK;
     EndIf
    Wend
   EndIf
  Wend
 EndIf

End Function

Function azimuth() ;-pa

 Local x#,y#,z#,x1#,y1#,z1#,zz#,theta1#,theta2#,ymin#,ymax#,theta0#;
 Local i,j;

 ymin = 2.0;
 ymax = -2.0;

 For j = 0 To Height-1
  If (debug And ((j Mod (Height/25)) = 0)) Write view
  For i = 0 To Width-1
   x = (2.0*i-Width)/Height/scale;
   y = (2.0*j-Height)/Height/scale;
   zz = x*x+y*y;
   z = 1.0-0.5*zz;
   If (z<-1.0)
    col(i,j) = BACK;
    If (doshade) shades(i,j) = 255;
   Else
    zz = Sqr(1.0-0.25*zz);
    x = x*zz;
    y = y*zz;
    x1 = clo*x+slo*sla*y+slo*cla*z;
    y1 = cla*y-sla*z;
    z1 = -slo*x+clo*sla*y+clo*cla*z;
    If (y1 < ymin) ymin = y1;
    If (y1 > ymax) ymax = y1;
    col(i,j) = planet0(x1,y1,z1);
    If (doshade) shades(i,j) = shade;
   EndIf
  Next
 Next

 If (hgrid <> 0.0) ;draw horisontal gridlines
  theta0 = 0.0
  While theta0>-90.0
   theta0=theta0-hgrid
   theta1 = theta0
   While theta1<90.0
    theta1=theta1+hgrid
    y = Sin(DEG2RAD*theta1*Rad2Deg) ;Sin(DEG2RAD*theta1)
    If (ymin <= y And y <= ymax)
     zz = Sqr(1-y*y);
     theta2=-Pi
     While theta2<Pi
      theta2=theta2+0.5/Width/scale
       x = Sin(theta2*Rad2Deg)*zz ;Sin(theta2)*zz
       z = Cos(theta2*Rad2Deg)*zz ;Cos(theta2)*zz
       x1 = clo*x-slo*z;
       y1 = slo*sla*x+cla*y+clo*sla*z;
       z1 = slo*cla*x-sla*y+clo*cla*z;
       If (z1<>-1.0)
        i = 0.5*(Height*scale*x1/Sqr(0.5+0.5*z1)+Width);
        j = 0.5*(Height*scale*y1/Sqr(0.5+0.5*z1)+Height);
        If (i>=0 And i<Width And j>=0 And j<Height) col(i,j) = BLACK;
       EndIf
      Wend
     EndIf
   Wend
  Wend
 EndIf

 If (vgrid <> 0.0) ;draw vertical gridlines
  theta2=-Pi
  While theta2<Pi
   theta2=theta2+0.5/Width/scale
   y = Sin(theta2*Rad2Deg) ;Sin(theta2)
   If (ymin <= y And y <= ymax)
    theta1 = 0.0
    While theta1<360.0
     theta1=theta1+vgrid
     x = Sin(DEG2RAD*theta1*Rad2Deg)*Cos(theta2*Rad2Deg);
     z = Cos(DEG2RAD*theta1*Rad2Deg)*Cos(theta2*Rad2Deg);
     ;x = Sin(DEG2RAD*theta1)*Cos(theta2);
     ;z = Cos(DEG2RAD*theta1)*Cos(theta2);
     x1 = clo*x-slo*z;
     y1 = slo*sla*x+cla*y+clo*sla*z;
     z1 = slo*cla*x-sla*y+clo*cla*z;
     If (z1<>-1.0)
      i = 0.5*(Height*scale*x1/Sqr(0.5+0.5*z1)+Width);
      j = 0.5*(Height*scale*y1/Sqr(0.5+0.5*z1)+Height);
      If (i>=0 And i<Width And j>=0 And j<Height) col(i,j) = BLACK;
     EndIf
    Wend
   EndIf
  Wend
 EndIf

End Function

Function conical() ;-pc

 Local k1#,c#,y2#,x#,y#,zz#,x1#,y1#,z1#
 Local theta1#,theta2#,ymin#,ymax#,cos2#,theta0#;
 Local i,j;

 ymin = 2.0;
 ymax = -2.0;

 If (lat>0)

  k1 = 1.0/Sin(lat*Rad2Deg) ;1.0/Sin(lat)
  c = k1*k1;
  y2 = Sqr(c*(1.0-Sin((lat/k1)*Rad2Deg))/(1.0+Sin((lat/k1)*Rad2Deg)));
  ;y2 = Sqr(c*(1.0-Sin(lat/k1))/(1.0+Sin(lat/k1)));

  For j = 0 To Height-1
   If (debug And ((j Mod (Height/25)) = 0)) Write view
   For i = 0 To Width-1
    x = (2.0*i-Width)/Height/scale;
    y = (2.0*j-Height)/Height/scale+y2;
    zz = x*x+y*y;
    If (zz=0.0) Then theta1 = 0.0 Else theta1 = k1*ATan2(x,y);
    If (theta1<-Pi Or theta1>Pi)
     col(i,j) = BACK;
     If (doshade) shades(i,j) = 255;
    Else
     theta1 = theta1 + longi-0.5*Pi ;theta1 is longitude
     theta2 = k1*ASin((zz-c)/(zz+c)) ;theta2 is latitude
     If (theta2 > 0.5*Pi Or theta2 < -0.5*Pi)
      col(i,j) = BACK;
      If (doshade) shades(i,j) = 255;
     Else
      cos2 = Cos(theta2*Rad2Deg) ;Cos(theta2)
      y = Sin(theta2*Rad2Deg) ;Sin(theta2)
      If (y < ymin) ymin = y;
      If (y > ymax) ymax = y;
      col(i,j) = planet0(Cos(theta1*Rad2Deg)*cos2,y,-Sin(theta1*Rad2Deg)*cos2);
      ;col(i,j) = planet0(Cos(theta1)*cos2,y,-Sin(theta1)*cos2);
      If (doshade) shades(i,j) = shade;
     EndIf
    EndIf
   Next
  Next

  If (hgrid <> 0.0) ;draw horisontal gridlines
   theta0 = 0.0
   While theta0>-90.0
    theta0=theta0-hgrid
    theta1 = theta0
    While theta1<90.0
    theta1=theta1+hgrid
     y = Sin(DEG2RAD*theta1*Rad2Deg) ;Sin(DEG2RAD*theta1)
     If (ymin <= y And y <= ymax)
      zz = Sqr(c*(1.0+Sin(DEG2RAD*(theta1/k1)*Rad2Deg))/(1.0-Sin(DEG2RAD*(theta1/k1)*Rad2Deg)));
      ;zz = Sqr(c*(1.0+Sin(DEG2RAD*theta1/k1))/(1.0-Sin(DEG2RAD*theta1/k1)));
      theta2=-Pi+longi
      While theta2<Pi+longi
       theta2=theta2+0.5/Width/scale
       z1 = theta2-longi;
       x1 = zz*Sin((z1/k1)*Rad2Deg) ;zz*Sin(z1/k1)
       y1 = zz*Cos((z1/k1)*Rad2Deg) ;zz*Cos(z1/k1)
       i = 0.5*(Height*scale*x1+Width);
       j = 0.5*(Height*scale*(y1-y2)+Height);
       If (i>=0 And i<Width And j>=0 And j<Height) col(i,j) = BLACK;
      Wend
     EndIf
    Wend
   Wend
  EndIf

  If (vgrid <> 0.0) ;draw vertical gridlines
   theta1=-0.5*Pi
   While theta1<0.5*Pi
    theta1=theta1+0.5/Width/scale
    y = Sin(theta1*Rad2Deg) ;Sin(theta1)
    If (ymin <= y And y <= ymax)
     zz = Sqr(c*(1.0+Sin(theta1/k1*Rad2Deg))/(1.0-Sin(theta1/k1*Rad2Deg)));
     ;zz = Sqr(c*(1.0+Sin(theta1/k1))/(1.0-Sin(theta1/k1)));
     theta0 = 0.0
     While theta0>-180.0+longi/DEG2RAD
      theta0=theta0-vgrid
      theta2 = theta0
      While theta2<180.0+longi/DEG2RAD
       theta2=theta2+vgrid
       z1 = DEG2RAD*theta2-longi
       x1 = zz*Sin(z1/k1*Rad2Deg) ;zz*Sin(z1/k1)
       y1 = zz*Cos(z1/k1*Rad2Deg) ;zz*Cos(z1/k1)
       i = 0.5*(Height*scale*x1+Width);
       j = 0.5*(Height*scale*(y1-y2)+Height);
       If (i>=0 And i<Width And j>=0 And j<Height) col(i,j) = BLACK;
      Wend
     Wend
    EndIf
   Wend
  EndIf

 Else

  k1 = 1.0/Sin(lat*Rad2Deg) ;1.0/Sin(lat)
  c = k1*k1;
  y2 = Sqr(c*(1.0-Sin(lat/k1*Rad2Deg))/(1.0+Sin(lat/k1*Rad2Deg)));
  ;y2 = Sqr(c*(1.0-Sin(lat/k1))/(1.0+Sin(lat/k1)));

  For j = 0 To Height-1
   If (debug And ((j Mod (Height/25)) = 0)) Write view
   For i = 0 To Width-1
    x = (2.0*i-Width)/Height/scale;
    y = (2.0*j-Height)/Height/scale-y2;
    zz = x*x+y*y;
    If (zz=0.0) Then theta1 = 0.0 Else theta1 = -k1*ATan2(x,-y);
    If (theta1<-Pi Or theta1>Pi)
     col(i,j) = 2;BACK;
     If (doshade) shades(i,j) = 255;
    Else
     theta1 = theta1 + longi-0.5*Pi ;theta1 is longitude
     theta2 = k1*ASin((zz-c)/(zz+c)) ;theta2 is latitude
     If (theta2 > 0.5*Pi Or theta2 < -0.5*Pi)
      col(i,j) = BACK;
      If (doshade) shades(i,j) = 255;
     Else
      cos2 = Cos(theta2*Rad2Deg) ;Cos(theta2)
      y = Sin(theta2*Rad2Deg) ;Sin(theta2)
      If (y < ymin) ymin = y;
      If (y > ymax) ymax = y;
      col(i,j) = planet0(Cos(theta1*Rad2Deg)*cos2,y,-Sin(theta1*Rad2Deg)*cos2);
      ;col(i,j) = planet0(Cos(theta1)*cos2,y,-Sin(theta1)*cos2);
      If (doshade) shades(i,j) = shade;
     EndIf
    EndIf
   Next
  Next

  If (hgrid <> 0.0) ;draw horisontal gridlines
   theta0 = 0.0
   While theta0>-90.0
    theta0=theta0-hgrid
    theta1 = theta0
    While theta1<90.0
     theta1=theta1+hgrid
     y = Sin(DEG2RAD*theta1*Rad2Deg) ;Sin(DEG2RAD*theta1)
     If (ymin <= y And y <= ymax)
      zz = Sqr(c*(1.0+Sin(DEG2RAD*theta1/k1*Rad2Deg))/(1.0-Sin(DEG2RAD*theta1/k1*Rad2Deg)));
      ;zz = Sqr(c*(1.0+Sin(DEG2RAD*theta1/k1))/(1.0-Sin(DEG2RAD*theta1/k1)));
      theta2=-Pi+longi
      While theta2<Pi+longi
       theta2=theta2+0.5/Width/scale
       z1 = theta2-longi;
       x1 = -zz*Sin(z1/k1*Rad2Deg) ;-zz*Sin(z1/k1)
       y1 = -zz*Cos(z1/k1*Rad2Deg) ;-zz*Cos(z1/k1)
       i = 0.5*(Height*scale*x1+Width);
       j = 0.5*(Height*scale*(y1+y2)+Height);
       If (i>=0 And i<Width And j>=0 And j<Height) col(i,j) = BLACK;
      Wend
     EndIf
    Wend
   Wend
  EndIf

  If (vgrid <> 0.0) ;draw vertical gridlines
   theta1=-0.5*Pi
   While theta1<0.5*Pi
    theta1=theta1+0.5/Width/scale
    y = Sin(theta1*Rad2Deg) ;Sin(theta1)
    If (ymin <= y And y <= ymax)
     zz = Sqr(c*(1.0+Sin(theta1/k1*Rad2Deg))/(1.0-Sin(theta1/k1*Rad2Deg)));
     ;zz = Sqr(c*(1.0+Sin(theta1/k1))/(1.0-Sin(theta1/k1)));
     theta0 = 0.0
     While theta0>-180.0+longi/DEG2RAD
      theta0=theta0-vgrid
      theta2 = theta0
      While theta2<180.0+longi/DEG2RAD
       theta2=theta2+vgrid
       z1 = DEG2RAD*theta2-longi;
       x1 = -zz*Sin(z1/k1*Rad2Deg) ;-zz*Sin(z1/k1)
       y1 = -zz*Cos(z1/k1*Rad2Deg) ;-zz*Cos(z1/k1)
       i = 0.5*(Height*scale*x1+Width);
       j = 0.5*(Height*scale*(y1+y2)+Height);
       If (i>=0 And i<Width And j>=0 And j<Height) col(i,j) = BLACK;
      Wend
     Wend
    EndIf
   Wend
  EndIf

 EndIf

End Function

Function heightfield()

 Local x#,y#,z#,x1#,y1#,z1#
 Local i,j;

 For j = 0 To Height-1
  If (debug And ((j Mod (Height/25)) = 0)) Write view
  For i = 0 To Width-1
   x = (2.0*i-Width)/Height/scale;
   y = (2.0*j-Height)/Height/scale;
   If (x*x+y*y>1.0)
    heights(i,j) = 0;
   Else
    z = Sqr(1.0-x*x-y*y);
    x1 = clo*x+slo*sla*y+slo*cla*z;
    y1 = cla*y-sla*z;
    z1 = -slo*x+clo*sla*y+clo*cla*z;
    heights(i,j) = 10000000*planet1(x1,y1,z1);
   EndIf
  Next
 Next

End Function

Function readmap()

 Local i,j;
 Local y#;
 Local c$,cline$ ;char

 Width = 47: Height = 21;

 For j = 0 To Height-1
  y = 0.5*7.5*(2.0*j-Height+1);
  y = Cos(DEG2RAD*y*Rad2Deg) ;Cos(DEG2RAD*y)
  weight(j) = Int(100.0*y+0.5);
 Next

 For j = 0 To Height-1 Step 2 ;j+=2
  Read cline ;read map data line
  For i = 0 To Width-1 Step 2 ;i+=2
   c = Mid(cline,1+(i/2),1) ;getchar()
   Select (c)
   Case ".": cl0(i,j) = -8;
   Case ",": cl0(i,j) = -4;
   Case ":": cl0(i,j) = -2;
   Case ";": cl0(i,j) = -1;
   Case "-": cl0(i,j) = 0;
   Case "*": cl0(i,j) = 1;
   Case "o": cl0(i,j) = 2;
   Case "O": cl0(i,j) = 4;
   Case "@": cl0(i,j) = 8;
   Default: Print "Wrong map symbol: "+c
   End Select
   If (i>0) cl0(i-1,j) = (cl0(i,j)+cl0(i-2,j))/2;
  Next
  If (debug) Print "line "+cline
 Next

 For j = 1 To Height-1 Step 2 ;j+=2
  For i = 0 To Width-1
   cl0(i,j) = (cl0(i,j-1)+cl0(i,j+1))/2;
  Next
 Next

End Function

Function search(file)

 Local y#,cos2#,theta1#,scale1#
 Local y2#,cos22#,theta12#;
 Local i,j,k,l,c,c1,c2,c3,errcount,errcount1,ascii$

 For j = 0 To Height-1
  y = 0.5*7.5*(2.0*j-Height+1);
  y = Sin(DEG2RAD*y*Rad2Deg) ;Sin(DEG2RAD*y)
  scale1 = Width/Height/Sqr(1.0-y*y)/Pi;
  cos2 = Sqr(1.0-y*y);
  y2 = 0.5*7.5*(2.0*j-Height+1.5);
  y2 = Sin(DEG2RAD*y2*Rad2Deg) ;Sin(DEG2RAD*y2)
  cos2 = Sqr(1.0-y2*y2);
  Depth = 3*(Int(log_2(scale1*Height)))+6;
  For i = 0 To Width-1
   theta1 = -0.5*Pi+Pi*(2.0*i-Width)/Width;
   theta12 = -0.5*Pi+Pi*(2.0*i+0.5-Width)/Width;
   c = 128+1000*planet1(Cos(theta1*Rad2Deg)*cos2,y,-Sin(theta1*Rad2Deg)*cos2);
   c1 = 128+1000*planet1(Cos(theta12*Rad2Deg)*cos2,y,-Sin(theta12*Rad2Deg)*cos2);
   c2 = 128+1000*planet1(Cos(theta1*Rad2Deg)*cos22,y2,-Sin(theta1*Rad2Deg)*cos22);
   c3 = 128+1000*planet1(Cos(theta12*Rad2Deg)*cos22,y2,-Sin(theta12*Rad2Deg)*cos22);
   ;c = 128+1000*planet1(Cos(theta1)*cos2,y,-Sin(theta1)*cos2);
   ;c1 = 128+1000*planet1(Cos(theta12)*cos2,y,-Sin(theta12)*cos2);
   ;c2 = 128+1000*planet1(Cos(theta1)*cos22,y2,-Sin(theta1)*cos22);
   ;c3 = 128+1000*planet1(Cos(theta12)*cos22,y2,-Sin(theta12)*cos22);
   c = (c+c1+c2+c3)/4.0;
   If (c<0) c = 0;
   If (c>255) c = 255;
   col(i,j) = c;
  Next
 Next

 For k=0 To Width-1
  For l=-20 To 20 Step 2 ;l+=2
   errcount = 0;

   For j = 0 To Height-1
    errcount1 = 0;
    For i = 0 To Width-1
     If (cl0(i,j)<0 And col((i+k) Mod Width,j) > 128-l)
      errcount1=errcount1-cl0(i,j);
     EndIf
     If (cl0(i,j)>0 And col((i+k) Mod Width,j) <= 128-l)
      errcount1=errcount1+cl0(i,j);
     EndIf
    Next
    errcount = errcount + weight(j)*errcount1;
   Next

   If (errcount < best)
    ascii="Errors: "+errcount+", parameters: -s "+rseed+" -l "+((360.0*k)/(Width+1))+" -i "+(M+l/1000.0)
    Print ascii
    best = errcount;
    For j = 0 To Height-1
     If j=0 ;first line of ascii map
      WriteByte file,13: WriteByte file,10 ;newline
      For i=1 To Len(ascii): WriteByte file,Asc(Mid(ascii,i,1)): Next
      WriteByte file,13: WriteByte file,10 ;newline
     EndIf
     For i = 0 To Width-1
      If (col((i+k) Mod Width,j) <= 128-l)
       WriteByte file,Asc(".") ;putchar('.');
      Else
       WriteByte file,Asc("O")
      EndIf
     Next
     WriteByte file,13: WriteByte file,10 ;newline
    Next
   EndIf
  Next
 Next

End Function

Function planet0(x#,y#,z#)
 ;x#,y#,z#

 Local alt#;
 Local colour;

 alt = planet1(x,y,z);

 If (altColors)

  Local snow# = 0.125;
  Local tree# = snow * 0.5;
  Local bare# = (tree + snow) / 2.0;
    
  If (latic)
   snow = snow - (0.13 * (y*y*y*y*y*y));
   bare = bare - (0.12 * (y*y*y*y*y*y));
   tree = tree - (0.11 * (y*y*y*y*y*y));
  EndIf

  If (alt > 0) ;Land
   If (alt > snow) ;Snow: White
    colour = WHITE;
   ElseIf (alt > bare) ;Snow: Grey - White
    colour = GREY0+Int((1+LAND1-GREY0) * (alt-bare)/(snow-bare));
    If (colour > LAND1) colour = LAND1;
   ElseIf (alt > tree) ;Bare: Brown - Grey
    colour = GREEN1+Int((1+GREY0-GREEN1) * (alt-tree)/(bare-tree));
    If (colour > GREY0) colour = GREY0;
   Else ;Green: Green - Brown
    colour = LAND0+Int((1+GREEN1-LAND0) * (alt)/(tree));
    If (colour > GREEN1) colour = GREEN1;
   EndIf
  Else ;Sea
   alt = alt/2;
   If (alt > snow) ;Snow: White
    colour = WHITE;
   ElseIf (alt > bare)
    colour = GREY0+Int((1+LAND1-GREY0) * (alt-bare)/(snow-bare));
    If (colour > LAND1) colour = LAND1;
   Else
    colour = BLUE1+Int((BLUE1-BLUE0+1)*(25*alt));
    If (colour<BLUE0) colour = BLUE0;
   EndIf
  EndIf

 Else ;calculate colour

  If (alt <=0.0) ;if below sea level then
   If (latic And y*y+alt >= 0.98)
    colour = LAND4 ;white if close to poles
   Else
    colour = BLUE1+Int((BLUE1-BLUE0+1)*(10*alt)) ;blue scale otherwise
    If (colour<BLUE0) colour = BLUE0;
   EndIf
  Else
   If (latic) alt = alt + 0.10204*y*y ;altitude adjusted with latitude
   If (alt >= 0.1) ;if high then
    colour = LAND4;
   Else ;else green to brown scale
    colour = LAND0+Int((LAND2-LAND0+1)*(10*alt));
    If (colour>LAND2) colour = LAND2;
   EndIf
  EndIf

 EndIf

 Return (colour);

End Function

Function planet#(a#,b#,c#,d#,as#,bs#,cs#,ds#,ax#,ay#,az#,bx#,by#,bz#,cx#,cy#,cz#,dx#,dy#,dz#,x#,y#,z#,level)
 ;a#,b#,c#,d# = altitudes of the 4 vertices
 ;as#,bs#,cs#,ds# = seeds of the 4 vertices
 ;ax#,ay#,az#,bx#,by#,bz#,cx#,cy#,cz#,dx#,dy#,dz# = vertex coordinates
 ;x#,y#,z# = goal point
 ;level = levels to go

 Local ssa#,ssb#,ssc#,ssd#,ssas#,ssbs#,sscs#,ssds#,ssax#,ssay#,ssaz#
 Local ssbx#,ssby#,ssbz#,sscx#,sscy#,sscz#,ssdx#,ssdy#,ssdz#;

 Local abx#,aby#,abz#,acx#,acy#,acz#,adx#,ady#,adz#;
 Local bcx#,bcy#,bcz#,bdx#,bdy#,bdz#,cdx#,cdy#,cdz#;
 Local lab#,lac#,lad#,lbc#,lbd#,lcd#;
 Local ex#,ey#,ez#,e#,es#,es1#,es2#,es3#;
 Local eax#,eay#,eaz#,epx#,epy#,epz#;
 Local ecx#,ecy#,ecz#,edx#,edy#,edz#;
 Local x1#,y1#,z1#,x2#,y2#,z2#,l1#,tmp#;

 If (level>0)

  If (level=11)
   ssa=a: ssb=b: ssc=c: ssd=d: ssas=as: ssbs=bs: sscs=cs: ssds=ds;
   ssax=ax: ssay=ay: ssaz=az: ssbx=bx: ssby=by: ssbz=bz;
   sscx=cx: sscy=cy: sscz=cz: ssdx=dx: ssdy=dy: ssdz=dz;
  EndIf

  abx = ax-bx: aby = ay-by: abz = az-bz;
  acx = ax-cx: acy = ay-cy: acz = az-cz;
  lab = abx*abx+aby*aby+abz*abz;
  lac = acx*acx+acy*acy+acz*acz;

  If (lab<lac)
   Return (planet(a,c,b,d,as,cs,bs,ds,ax,ay,az,cx,cy,cz,bx,by,bz,dx,dy,dz,x,y,z,level));
  Else
   adx = ax-dx: ady = ay-dy: adz = az-dz;
   lad = adx*adx+ady*ady+adz*adz;
   If (lab<lad)
    Return (planet(a,d,b,c,as,ds,bs,cs,ax,ay,az,dx,dy,dz,bx,by,bz,cx,cy,cz,x,y,z,level));
   Else
    bcx = bx-cx: bcy = by-cy: bcz = bz-cz;
    lbc = bcx*bcx+bcy*bcy+bcz*bcz;
    If (lab<lbc)
     Return (planet(b,c,a,d,bs,cs,as,ds,bx,by,bz,cx,cy,cz,ax,ay,az,dx,dy,dz,x,y,z,level));
    Else
     bdx = bx-dx: bdy = by-dy: bdz = bz-dz;
     lbd = bdx*bdx+bdy*bdy+bdz*bdz;
     If (lab<lbd)
      Return (planet(b,d,a,c,bs,ds,as,cs,bx,by,bz,dx,dy,dz,ax,ay,az,cx,cy,cz,x,y,z,level));
     Else
      cdx = cx-dx: cdy = cy-dy: cdz = cz-dz;
      lcd = cdx*cdx+cdy*cdy+cdz*cdz;
      If (lab<lcd)
       Return (planet(c,d,a,b,cs,ds,as,bs,cx,cy,cz,dx,dy,dz,ax,ay,az,bx,by,bz,x,y,z,level));
      Else
       es = rand2(as,bs);
       es1 = rand2(es,es);
       es2 = 0.5+0.1*rand2(es1,es1);
       es3 = 1.0-es2;
       If (ax=bx) ;very unlikely to ever happen
        ex = 0.5*ax+0.5*bx: ey = 0.5*ay+0.5*by: ez = 0.5*az+0.5*bz;
       ElseIf (ax<bx)
        ex = es2*ax+es3*bx: ey = es2*ay+es3*by: ez = es2*az+es3*bz;
       Else
        ex = es3*ax+es2*bx: ey = es3*ay+es2*by: ez = es3*az+es2*bz;
       EndIf
       If (lab>1.0) lab = lab^0.75 ;pow(lab,0.75);
       e = 0.5*(a+b)+es*dd1*Abs(a-b)+es1*dd2*(lab^POW) ;pow(lab,POW);
       eax = ax-ex: eay = ay-ey: eaz = az-ez;
       epx =  x-ex: epy =  y-ey: epz =  z-ez;
       ecx = cx-ex: ecy = cy-ey: ecz = cz-ez;
       edx = dx-ex: edy = dy-ey: edz = dz-ez;
       If ((eax*ecy*edz+eay*ecz*edx+eaz*ecx*edy-eaz*ecy*edx-eay*ecx*edz-eax*ecz*edy)*(epx*ecy*edz+epy*ecz*edx+epz*ecx*edy-epz*ecy*edx-epy*ecx*edz-epx*ecz*edy)>0.0)
        Return (planet(c,d,a,e,cs,ds,as,es,cx,cy,cz,dx,dy,dz,ax,ay,az,ex,ey,ez,x,y,z,level-1));
       Else
        Return (planet(c,d,b,e,cs,ds,bs,es,cx,cy,cz,dx,dy,dz,bx,by,bz,ex,ey,ez,x,y,z,level-1));
       EndIf
      EndIf
     EndIf
    EndIf
   EndIf 
  EndIf

 Else

  If (doshade)
   x1 = 0.25*(ax+bx+cx+dx);
   x1 = a*(x1-ax)+b*(x1-bx)+c*(x1-cx)+d*(x1-dx);
   y1 = 0.25*(ay+by+cy+dy);
   y1 = a*(y1-ay)+b*(y1-by)+c*(y1-cy)+d*(y1-dy);
   z1 = 0.25*(az+bz+cz+dz);
   z1 = a*(z1-az)+b*(z1-bz)+c*(z1-cz)+d*(z1-dz);
   l1 = Sqr(x1*x1+y1*y1+z1*z1);
   If (l1=0.0) l1 = 1.0;
   tmp = Sqr(1.0-y*y);
   If (tmp<0.0001) tmp = 0.0001;
   x2 = x*x1+y*y1+z*z1;
   y2 = -x*y/tmp*x1+tmp*y1-z*y/tmp*z1;
   z2 = -z/tmp*x1+x/tmp*z1;
   shade = Int((-Sin(Pi*shade_angle/180.0)*y2-Cos(Pi*shade_angle/180.0)*z2)/l1*48.0+128.0);
   If (shade<10) shade = 10;
   If (shade>255) shade = 255;
  EndIf
  Return ((a+b+c+d)/4);

 EndIf

End Function

Function planet1#(x#,y#,z#)
 ;x#,y#,z#

 Local abx#,aby#,abz#,acx#,acy#,acz#,adx#,ady#,adz#,apx#,apy#,apz#;
 Local bax#,bay#,baz#,bcx#,bcy#,bcz#,bdx#,bdy#,bdz#,bpx#,bpy#,bpz#;

 abx = ssbx-ssax: aby = ssby-ssay: abz = ssbz-ssaz;
 acx = sscx-ssax: acy = sscy-ssay: acz = sscz-ssaz;
 adx = ssdx-ssax: ady = ssdy-ssay: adz = ssdz-ssaz;
 apx = x-ssax: apy = y-ssay: apz = z-ssaz;

 If ((adx*aby*acz+ady*abz*acx+adz*abx*acy-adz*aby*acx-ady*abx*acz-adx*abz*acy)*(apx*aby*acz+apy*abz*acx+apz*abx*acy-apz*aby*acx-apy*abx*acz-apx*abz*acy)>0.0);
  ;p is on same side of abc as d
  If ((acx*aby*adz+acy*abz*adx+acz*abx*ady-acz*aby*adx-acy*abx*adz-acx*abz*ady)*(apx*aby*adz+apy*abz*adx+apz*abx*ady-apz*aby*adx-apy*abx*adz-apx*abz*ady)>0.0);
   ;p is on same side of abd as c
   If ((abx*ady*acz+aby*adz*acx+abz*adx*acy-abz*ady*acx-aby*adx*acz-abx*adz*acy)*(apx*ady*acz+apy*adz*acx+apz*adx*acy-apz*ady*acx-apy*adx*acz-apx*adz*acy)>0.0);
    ;p is on same side of acd as b
    bax = -abx: bay = -aby: baz = -abz;
    bcx = sscx-ssbx: bcy = sscy-ssby: bcz = sscz-ssbz;
    bdx = ssdx-ssbx: bdy = ssdy-ssby: bdz = ssdz-ssbz;
    bpx = x-ssbx: bpy = y-ssby: bpz = z-ssbz;
    If ((bax*bcy*bdz+bay*bcz*bdx+baz*bcx*bdy-baz*bcy*bdx-bay*bcx*bdz-bax*bcz*bdy)*(bpx*bcy*bdz+bpy*bcz*bdx+bpz*bcx*bdy-bpz*bcy*bdx-bpy*bcx*bdz-bpx*bcz*bdy)>0.0);
     ;p is on same side of bcd as a
     ;Hence, p is inside tetrahedron
     Return (planet(ssa,ssb,ssc,ssd,ssas,ssbs,sscs,ssds,ssax,ssay,ssaz,ssbx,ssby,ssbz,sscx,sscy,sscz,ssdx,ssdy,ssdz,x,y,z,11));
    EndIf
   EndIf
  EndIf
 EndIf

 ;otherwise
 Return (planet(M,M,M,M,r1,r2,r3,r4,0.0,0.0,3.01,0.0,Sqr(8.0)+0.01*r1*r1,-1.02+0.01*r2*r3,-Sqr(6.0)-0.01*r3*r3,-Sqr(2.0)-0.01*r4*r4,-1.02+0.01*r1*r2,Sqr(6.0)-0.01*r2*r2,-Sqr(2.0)-0.01*r3*r3,-1.02+0.01*r1*r3,x,y,z,Depth));
 ;M,M,M,M = initial altitude is M on all corners of tetrahedron
 ;r1,r2,r3,r4 = same seed set is used in every call
 ;coordinates of vertices
 ;x,y,z = coordinates of point we want colour of
 ;Depth = subdivision depth

End Function

Function rand2#(p#,q#)
 ;p#,q#
 ;random number generator taking two seeds
 ;rand2(p,q) = rand2(q,p) is important

 Local r#;
 r = (p+3.14159265)*(q+3.14159265);
 Return (2.0*(r-Int(r))-1.0);

End Function

Function printppm(filename$)
 ;prints picture in PPM (portable pixel map) format

 Local i,j,c,s,file,ascii$;

 file=WriteFile(filename)
 If Not file Return False ;write fail

 ascii="P6"+Chr(10) ;"P6\n"
 For i=1 To Len(ascii): WriteByte file,Asc(Mid(ascii,i,1)): Next
 ascii="#fractal planet image"+Chr(10) ;"#title\n"
 For i=1 To Len(ascii): WriteByte file,Asc(Mid(ascii,i,1)): Next
 ascii=Str(Width)+" "+Str(Height)+" "+Str(255)+Chr(10) ;"info ncols\n"
 For i=1 To Len(ascii): WriteByte file,Asc(Mid(ascii,i,1)): Next

 If (doshade)
  For j=0 To Height-1
   For i=0 To Width-1
    s = shades(i,j);
    c = s*rtable(col(i,j))/150;
    If (c>255) c=255;
    WriteByte file,c;
    c = s*gtable(col(i,j))/150;
    If (c>255) c=255;
    WriteByte file,c;
    c = s*btable(col(i,j))/150;
    If (c>255) c=255;
    WriteByte file,c;
   Next
  Next
 Else
  For j=0 To Height-1
   For i=0 To Width-1
    WriteByte file,rtable(col(i,j));
    WriteByte file,gtable(col(i,j));
    WriteByte file,btable(col(i,j));
   Next
  Next
 EndIf

 CloseFile file;

End Function

Function printppmBW(filename$)
 ;prints picture in b/w PPM format

 Local i,j,c,file,ascii$;

 file=WriteFile(filename)
 If Not file Return False ;write fail

 ascii="P6"+Chr(10) ;"P6\n"
 For i=1 To Len(ascii): WriteByte file,Asc(Mid(ascii,i,1)): Next
 ascii="#fractal planet image"+Chr(10) ;"#title\n"
 For i=1 To Len(ascii): WriteByte file,Asc(Mid(ascii,i,1)): Next
 ascii=Str(Width)+" "+Str(Height)+" "+Str(2)+Chr(10) ;"info ncols\n"
 For i=1 To Len(ascii): WriteByte file,Asc(Mid(ascii,i,1)): Next

  For j=0 To Height-1
   For i=0 To Width-1
   If (col(i,j) = BLACK) Then c=0 Else c=1;
   WriteByte file,c;
   WriteByte file,c;
   WriteByte file,c;
  Next
 Next

 CloseFile file;

End Function

Function printbmp(filename$)
 ;prints picture in BMP format

 Local i,j,k,c,s,W1,file;

 file=WriteFile(filename)
 If Not file Return False ;write fail

 W1 = (3*Width+3);
 W1 = W1 - W1 Mod 4;
 s = 54+W1*Height ;file size

 WriteShort file,(Asc("M") Shl 8)+Asc("B") ;file type "BM"
 WriteInt file,s ;bfSize
 WriteShort file,0 ;bfReserved1
 WriteShort file,0 ;bfReserved2
 WriteInt file,54 ;bfOffBits, offset to data

 WriteInt file,40 ;size of infoheader
 WriteInt file,Width ;biWidth
 WriteInt file,Height ;biHeight
 WriteShort file,1 ;no. of planes = 1
 WriteShort file,24 ;bpp
 WriteInt file,0 ;no compression
 WriteInt file,0 ;image size (unspecified)
 WriteInt file,32 Shl 8 ;h. pixels/m
 WriteInt file,32 Shl 8 ;v. pixels/m
 WriteInt file,0 ;colours used (unspecified)
 WriteInt file,0 ;important colours (all)

 If (doshade)
  For k=0 To Height-1 ;j--
   j=Height-1-k
   For i=0 To Width-1
    s = shades(i,j);
    c = s*btable(col(i,j))/150;
    If (c>255) c=255;
    WriteByte file,c;
    c = s*gtable(col(i,j))/150;
    If (c>255) c=255;
    WriteByte file,c
    c = s*rtable(col(i,j))/150;
    If (c>255) c=255;
    WriteByte file,c
   Next
   For i=3*Width To W1-1: WriteByte file,0: Next
  Next
 Else
  For k=0 To Height-1
   j=Height-1-k
   For i=0 To Width-1
    WriteByte file,btable(col(i,j));
    WriteByte file,gtable(col(i,j));
    WriteByte file,rtable(col(i,j));
   Next
   For i=3*Width To W1-1: WriteByte file,0: Next
  Next
 EndIf

 CloseFile file;

End Function

Function printbmpBW(filename$)
 ;prints picture in b/w BMP format

 Local i,j,k,c,s,W1;

 file=WriteFile(filename)
 If Not file Return False ;write fail

 W1 = (Width+31);
 W1 = W1 - W1 Mod 32;
 s = 62+(W1*Height)/8 ;file size

 WriteShort file,(Asc("M") Shl 8)+Asc("B") ;file type "BM"
 WriteInt file,s ;bfSize
 WriteShort file,0 ;bfReserved1
 WriteShort file,0 ;bfReserved2
 WriteInt file,62 ;bfOffBits, offset to data

 WriteInt file,40 ;size of infoheader
 WriteInt file,Width ;biWidth
 WriteInt file,Height ;biHeight
 WriteShort file,1 ;no. of planes = 1
 WriteShort file,1 ;bpp
 WriteInt file,0 ;no compression
 WriteInt file,0 ;image size (unspecified)
 WriteInt file,32 Shl 8 ;h. pixels/m
 WriteInt file,32 Shl 8 ;v. pixels/m
 WriteInt file,2 ;colours used (unspecified)
 WriteInt file,2 ;important colours (all)

 WriteInt file,0 ;colour 0 = black
 WriteInt file,$FFFFFFFF ;colour 1 = white

 For k=0 To Height-1 ;j--
  j=Height-1-k
  For i=0 To W1-1 Step 8 ;i+=8
   If (i<Width And col(i,j) <> BLACK) c=128 Else c=0;
   If (i+1<Width And col(i+1,j) <> BLACK) c=c+64;
   If (i+2<Width And col(i+2,j) <> BLACK) c=c+32;
   If (i+3<Width And col(i+3,j) <> BLACK) c=c+16;
   If (i+4<Width And col(i+4,j) <> BLACK) c=c+8;
   If (i+5<Width And col(i+5,j) <> BLACK) c=c+4;
   If (i+6<Width And col(i+6,j) <> BLACK) c=c+2;
   If (i+7<Width And col(i+7,j) <> BLACK) c=c+1;
   WriteByte file,c;
  Next
 Next

 CloseFile file;

End Function

Function printheights(filename$)

 Local i,j,k,c,s,W1,file;

 file=WriteFile(filename)
 If Not file Return False ;write fail

 W1 = (3*Width+3);
 W1 = W1 - W1 Mod 4;
 s = 54+W1*Height ;file size

 WriteShort file,(Asc("M") Shl 8)+Asc("B") ;file type "BM"
 WriteInt file,s ;bfSize
 WriteShort file,0 ;bfReserved1
 WriteShort file,0 ;bfReserved2
 WriteInt file,54 ;bfOffBits, offset to data

 WriteInt file,40 ;size of infoheader
 WriteInt file,Width ;biWidth
 WriteInt file,Height ;biHeight
 WriteShort file,1 ;no. of planes = 1
 WriteShort file,24 ;bpp
 WriteInt file,0 ;no compression
 WriteInt file,0 ;image size (unspecified)
 WriteInt file,32 Shl 8 ;h. pixels/m
 WriteInt file,32 Shl 8 ;v. pixels/m
 WriteInt file,0 ;colours used (unspecified)
 WriteInt file,0 ;important colours (all)

 For j=0 To Height-1
  For i=0 To Width-1
   k=Height-1-j ;fprintf(filename,"%d ",heights[i][j])
   WriteByte file,(heights(i,k) And $0000FF);
   WriteByte file,((heights(i,k) And $00FF00) Shr 8);
   WriteByte file,(heights(i,k) Shr 16);
  Next
  For i=3*Width To W1-1: WriteByte file,0: Next
 Next

 CloseFile file;

End Function

Function log_2#(x#)
 Return (Log(x)/Log(2.0))
End Function

Function print_error()

 Print "Usage: [options]";
 Print " -?               (or any illegal option) Output this text";
 Print " -s seed          Seed as number between 0.0 and 1.0";
 Print " -w width         Width in pixels, default = 640";
 Print " -h height        Height in pixels, default = 480";
 Print " -m magnification Magnification, default = 1.0";
 Print " -o output_file   Output file, default = planet.bmp";
 Print " -l longitude     Longitude of centre in degrees, default = 0.0";
 Print " -L latitude      Latitude of centre in degrees, default = 0.0";
 Print " -g gridsize      Vertical gridsize in degrees, default = 0.0 (no grid)";
 Print " -G gridsize      Horisontal gridsize in degrees, default = 0.0 (no grid)";
 Print " -i init_alt      Initial altitude, default = 0.3";
 Print " -c               Colour depends on latitude, default = only altitude";
 Print " -C               Use lighter colours (original scheme only)";
 Print " -N nocols        Use nocols number of colours, default = 256, minimum = 5";
 Print " -a               Switch to alternative colour scheme";
 Print " -M file          Read colour definitions from file";
 Print " -O               Produce a black and white outline map";
 Print " -E               Trace the edges of land in black on colour map";
 Print " -B               Use 'bumpmap' shading";
 Print " -A angle         Angle of 'light' in bumpmap shading";
 Print " -P               Use PPM file format, default = BMP";
 Print " -b               Reverse background colour, default = black";
 Print " -V number        Distance contribution to variation, default = 0.03";
 Print " -v number        Altitude contribution to variation, default = 0.4";
 Print " -pprojection     Projection:";
 Print "  m = Mercator (default), p = Peters, q = Square, s = Stereographic";
 Print "  o = Orthographic, g = Gnomonic, a = Area preserving azimuthal";
 Print "  c = Conical (conformal), M = Mollweide, S = Sinusoidal";
 Print "  h = Heightfield, f = Find match, see manual";

 ;With the -pf option a map must be given on standard input.
 ;This map is 11 lines of 24 characters. The characters are:
 ; . : very strong preference for water (value=8)
 ; , : strong preference for water (value=4)
 ; : : preference for water (value=2)
 ; ; : weak preference for water (value=1)
 ; - : don't care (value=0)
 ; * : weak preference for land (value=1)
 ; o : preference for land (value=2)
 ; O : strong preference for land (value=4)
 ; @ : very strong preference for land (value=8)
 ;
 ;Each point on the map corresponds to a point on a 15� grid.
 ;
 ;The program tries seeds starting from the specified and
 ;successively outputs the seed (and rotation) of the best
 ;current match, together with a small map of this.
 ;This is all ascii, no bitmap is produced.

End Function

Comments

2006

It would not open your mars part and work.

markcw

2006

Hi NA, it is designed to work without the mars.rgb color file. You have to wait a while before you see anything. :)

I have updated the code now to show the -B option (bumpmapping).

Also i have produced a cut-down version specifically to map onto a 3d sphere in blitz3d using the squarep projection. It keeps all the main functional code so you can still play around with the important stuff.

edit: updated this, now it's simpler and bumpmap is smoother.
edit2: last update. added new color schemes and changed the bumpmap level a bit (it was too smooth).
edit3: ok, i updated this again because i noticed the top row of pixels was always black if graphics3d (windowed) was called before the map was saved. very strange.

;planet generating program in 3d, using squarep

;these values can be changed to change world characteristics

Global debug = 1 ;show debug info
Global colorscheme = 1 ;specify colour scheme 0..5, default = 0
Global altColors = 0 ;Switch to alternative colour scheme

Global doshade = 0 ;Use 'bumpmap' shading
Global shade_angle# = 0 ;Angle of 'light' in bumpmap shading, eg. 150.0

Global Width = 512 ;Width in pixels, default = 256
Global Height = 256 ;Height in pixels, default = 128

Global rseed# = 0.2 ;Seed as number between 0.0..1.0, eg. 0.3
Global M# = 0.19 ;Initial altitude, default = 0.3
Global dd1# = 0.15 ;Weight for altitude difference, default = 0.4
Global dd2# = 0.03 ;Weight for distance, default = 0.03
Global POW# = 0.2 ;Power for distance function, default = 0.47

Global scale# = 1.0 ;Magnification, default = 1.0
Global longi# = 0.0 ;Longitude of centre in degrees, default = 0.0
Global lat# = 0.0 ;Latitude of centre in degrees, default = 0.0
Global vgrid# = 0.0 ;Vertical gridsize in degrees, 0.0 = no grid
Global hgrid# = 0.0 ;Horisontal gridsize in degrees, 0.0 = no grid

Global latic = 0 ;Colour depends on latitude, default = only altitude
Global lighter = 0 ;Use lighter colours (original scheme only)

;globals used in calculations
Global BLACK=0
Global WHITE=1
Global BLUE0=2
Global BLUE1, LAND0, LAND1, LAND2, LAND4;
Global GREEN1, BROWN0, GREY0;
Global BACK = WHITE ;BLACK

Global shade;
Global Depth ;depth of subdivisions
Global r1#,r2#,r3#,r4# ;seeds
Global DEG2RAD#=0.0174532918661 ;Pi/180
Global Rad2Deg#=(180/Float(Pi))

;arrays used in calculations
Dim col(256,128)
Dim shades(256,128)
Dim rtable(256)
Dim gtable(256)
Dim btable(256)
Dim colors(9,3)

;set colors array from data
Select colorscheme
 Case 1: Restore atlas_colors
 Case 2: Restore mars_colors
 Case 3: Restore grey_colors
 Case 4: Restore moon_colors
 Case 5: Restore cyan_colors
 Default: Restore colors
End Select
For i=0 To 8 ;copy the 9 base colors
 For j=0 To 2: Read colors(i,j): Next
Next

Dim col(Width,Height) ;map colors
If (doshade) Dim shades(Width,Height) ;map shading

If (debug)
 Print "> -s "+rseed+" -i "+M+" -v "+dd1+" -V "+dd2+" -P "+POW
 Print "> -w "+Width+" -h "+Height+" -B "+doshade
EndIf

;init calculations
If (longi>180) longi = longi - 360;
longi = longi*DEG2RAD;
lat = lat*DEG2RAD;
Depth = 3*(Int(log_2(scale*Height)))+6;
r1 = rseed;
r1 = rand2(r1,r1);
r2 = rand2(r1,r1);
r3 = rand2(r1,r2);
r4 = rand2(r2,r3);

setcolours() ;set rgb arrays

If (debug) Print "+----+----+----+----+----+"

squarep() ;Square projection (equidistant latitudes)

If (doshade) smoothshades() ;bumpmap shading

;draw map
image=CreateImage(Width,Height)
LockBuffer ImageBuffer(image)
For y=0 To ImageHeight(image)-1
 For x=0 To ImageWidth(image)-1
  r=rtable(col(x,y))
  g=gtable(col(x,y))
  b=btable(col(x,y))
  If doshade ;bumpmap
   s = (shades(x,y)/3)+85 ;smooth bumpmap
   r = s*rtable(col(x,y))/150;
   If (r>255) r=255;
   g = s*gtable(col(x,y))/150;
   If (g>255) g=255;
   b = s*btable(col(x,y))/150;
   If (b>255) b=255;
  EndIf
  rgb=(r Shl 16)+(g Shl 8)+b
  WritePixelFast x,y,rgb,ImageBuffer(image)
 Next
Next
UnlockBuffer ImageBuffer(image)
SetBuffer BackBuffer()

SaveImage(image,"planet.bmp") ;save the result

;init 3d world
Graphics3D 640,480,0,2
SetBuffer BackBuffer()

camera=CreateCamera()
PositionEntity camera,0,0,-2
light=CreateLight()
RotateEntity light,45,45,0

ball=CreateSphere(20)
texture=LoadTexture("planet.bmp") ;load map
EntityTexture ball,texture

;main loop
While Not KeyHit(1) Or KeyHit(57) ;esc or space keys
 RenderWorld

 TurnEntity ball,0,0.1,0.07

 If KeyHit(17) wf=Not wf : WireFrame wf ;W key

 Flip
Wend
End

;data labels

.colors ;colorscheme 0
Data 0,0,255 ;Dark blue depths
Data 0,128,255 ;Light blue shores
Data 0,255,0 ;Light green lowlands
Data 64,192,16 ;Dark green highlands
Data 64,192,16 ;Dark green Mountains
Data 128,128,32 ;Brown stoney peaks
Data 255,255,255 ;White - peaks
Data 0,0,0 ;Black - Space
Data 0,0,0 ;Black - Lines

.atlas_colors ;colorscheme 1
Data 0,0,192 ;Dark blue depths
Data 0,128,255 ;Light blue shores
Data 0,96,0 ;Dark green Lowlands
Data 0,224,0 ;Light green Highlands
Data 128,176,0 ;Brown mountainsides
Data 128,128,128 ;Grey stoney peaks
Data 255,255,255 ;White - peaks
Data 0,0,0 ;Black - Space
Data 0,0,0 ;Black - Lines

.mars_colors ;colorscheme 2
Data 55,28,75 ;0
Data 136,92,80 ;1
Data 148,100,82 ;2
Data 232,164,90 ;3
Data 248,180,92 ;4
Data 165,128,108 ;5
Data 255,255,255 ;6
Data 0,0,0 ;7
Data 0,0,0 ;8

.grey_colors ;colorscheme 3
Data 96,96,96 ;0
Data 152,152,152 ;2
Data 144,144,144 ;1
Data 168,168,168 ;3
Data 160,160,160 ;4
Data 144,144,144 ;5
Data 128,128,128 ;6
Data 0,0,0 ;7
Data 0,0,0 ;8

.moon_colors ;colorscheme 4
Data 112,112,98 ;0
Data 168,168,147 ;2
Data 160,160,140 ;1
Data 184,184,161 ;3
Data 176,176,154 ;4
Data 160,160,140 ;5
Data 144,144,126 ;6
Data 0,0,0 ;7
Data 0,0,0 ;8

.cyan_colors ;colorscheme 5
Data 74,112,150 ;0
Data 126,168,210 ;2
Data 120,160,200 ;1
Data 138,184,230 ;3
Data 132,176,220 ;4
Data 120,160,200 ;5
Data 108,144,180 ;6
Data 0,0,0 ;7
Data 0,0,0 ;8

;functions

Function min(x,y)
 If x<y Then Return x Else Return y
End Function

Function max(x,y)
 If x<y Then Return y Else Return x
End Function

Function setcolours()

 Local i;
 Local crow;
 Local nocols = 256

 If (debug)
  For crow = 0 To 8
   Write "c"+crow+"="+colors(crow,0)+","+colors(crow,1)+","+colors(crow,2)+" ";
   If Not (crow+1) Mod 3 Print " "
  Next
 EndIf

 If (altColors)

  If (nocols < 8) nocols = 8;
    
  LAND0 = max(nocols / 4, BLUE0 + 1);
  BLUE1 = LAND0 - 1;
  GREY0 = nocols - (nocols / 8);
  GREEN1 = min(LAND0 + (nocols / 2), GREY0 - 2);
  BROWN0 = (GREEN1 + GREY0) / 2;
  LAND1 = nocols - 1;

  rtable(BLACK) = colors(7,0);
  gtable(BLACK) = colors(7,0);
  btable(BLACK) = colors(7,0);
  rtable(WHITE) = colors(6,0);
  gtable(WHITE) = colors(6,1);
  btable(WHITE) = colors(6,2);
  rtable(BLUE0) = colors(0,0);
  gtable(BLUE0) = colors(0,1);
  btable(BLUE0) = colors(0,2);

  For i=BLUE0+1 To BLUE1;
   rtable(i) = (colors(0,0)*(BLUE1-i)+colors(1,0)*(i-BLUE0))/(BLUE1-BLUE0);
   gtable(i) = (colors(0,1)*(BLUE1-i)+colors(1,1)*(i-BLUE0))/(BLUE1-BLUE0);
   btable(i) = (colors(0,2)*(BLUE1-i)+colors(1,2)*(i-BLUE0))/(BLUE1-BLUE0);
  Next
  For i=LAND0 To GREEN1-1
   rtable(i) = (colors(2,0)*(GREEN1-i)+colors(3,0)*(i-LAND0))/(GREEN1-LAND0);
   gtable(i) = (colors(2,1)*(GREEN1-i)+colors(3,1)*(i-LAND0))/(GREEN1-LAND0);
   btable(i) = (colors(2,2)*(GREEN1-i)+colors(3,2)*(i-LAND0))/(GREEN1-LAND0);
  Next
  For i=GREEN1 To BROWN0-1
   rtable(i) = (colors(3,0)*(BROWN0-i)+colors(4,0)*(i-GREEN1))/(BROWN0-GREEN1);
   gtable(i) = (colors(3,1)*(BROWN0-i)+colors(4,1)*(i-GREEN1))/(BROWN0-GREEN1);
   btable(i) = (colors(3,2)*(BROWN0-i)+colors(4,2)*(i-GREEN1))/(BROWN0-GREEN1);
  Next
  For i=BROWN0 To GREY0-1
   rtable(i) = (colors(4,0)*(GREY0-i)+colors(5,0)*(i-BROWN0))/(GREY0-BROWN0);
   gtable(i) = (colors(4,1)*(GREY0-i)+colors(5,1)*(i-BROWN0))/(GREY0-BROWN0);
   btable(i) = (colors(4,2)*(GREY0-i)+colors(5,2)*(i-BROWN0))/(GREY0-BROWN0);
  Next
  For i=GREY0 To nocols-1
   rtable(i) = (colors(5,0)*(nocols-i)+(colors(6,0)+1)*(i-GREY0))/(nocols-GREY0);
   gtable(i) = (colors(5,1)*(nocols-i)+(colors(6,1)+1)*(i-GREY0))/(nocols-GREY0);
   btable(i) = (colors(5,2)*(nocols-i)+(colors(6,2)+1)*(i-GREY0))/(nocols-GREY0);
  Next

 Else

  rtable(BLACK) = 0;
  gtable(BLACK) = 0;
  btable(BLACK) = 0;
  rtable(WHITE) = 255;
  gtable(WHITE) = 255;
  btable(WHITE) = 255;

  Local r, c;
  Local x#;

  While (lighter>0)
   For r = 0 To 7-1
    For c = 0 To 3-1
     x = Sqr(Float(colors(r,c))/256.0);
     colors(r,c) = Int(240.0*x+16);
    Next
   Next
   lighter=lighter-1
  Wend

  BLUE1 = (nocols-4)/2+BLUE0;

  If (BLUE1=BLUE0)
   rtable(BLUE0) = colors(0,0);
   gtable(BLUE0) = colors(0,1);
   btable(BLUE0) = colors(0,2);
  Else
   For i=BLUE0 To BLUE1;
    rtable(i) = (colors(0,0)*(BLUE1-i)+colors(1,0)*(i-BLUE0))/(BLUE1-BLUE0);
    gtable(i) = (colors(0,1)*(BLUE1-i)+colors(1,1)*(i-BLUE0))/(BLUE1-BLUE0);
    btable(i) = (colors(0,2)*(BLUE1-i)+colors(1,2)*(i-BLUE0))/(BLUE1-BLUE0);
   Next
  EndIf
  LAND0 = BLUE1+1: LAND2 = nocols-2: LAND1 = (LAND0+LAND2+1)/2;
  For i=LAND0 To LAND1-1
   rtable(i) = (colors(2,0)*(LAND1-i)+colors(3,0)*(i-LAND0))/(LAND1-LAND0);
   gtable(i) = (colors(2,1)*(LAND1-i)+colors(3,1)*(i-LAND0))/(LAND1-LAND0);
   btable(i) = (colors(2,2)*(LAND1-i)+colors(3,2)*(i-LAND0))/(LAND1-LAND0);
  Next

  If (LAND1=LAND2)
   rtable(LAND1) = colors(4,0);
   gtable(LAND1) = colors(4,1);
   btable(LAND1) = colors(4,2);
  Else
   For i=LAND1 To LAND2;
    rtable(i) = (colors(4,0)*(LAND2-i)+colors(5,0)*(i-LAND1))/(LAND2-LAND1);
    gtable(i) = (colors(4,1)*(LAND2-i)+colors(5,1)*(i-LAND1))/(LAND2-LAND1);
    btable(i) = (colors(4,2)*(LAND2-i)+colors(5,2)*(i-LAND1))/(LAND2-LAND1);
   Next
  EndIf
  LAND4 = nocols-1;
  rtable(LAND4) = colors(6,0);
  gtable(LAND4) = colors(6,1);
  btable(LAND4) = colors(6,2);

 EndIf

End Function

Function smoothshades()

 Local i,j;

 For i=0 To Width-3 ;i<Width-2
  For j=0 To Height-3;
   shades(i,j) = (4*shades(i,j)+2*shades(i,j+1)+2*shades(i+1,j)+shades(i+1,j+2)+4)/9;
  Next
 Next

End Function

Function squarep() ;-pq

 Local y#,scale1#,theta1#,cos2#,theta0#;
 Local k,i,j;

 k = Int(lat*Width*scale/Pi);

 For j = 0 To Height-1
  If (debug And ((j Mod (Height/25)) = 0)) Write "q"
  y = (2.0*(j-k)-Height)/Width/scale*Pi;
  If (Abs(y)>0.5*Pi) ;(Abs(y)>=0.5*Pi)
   For i = 0 To Width-1
    col(i,j) = BACK;
    If (doshade) shades(i,j) = 255;
   Next
  Else
   cos2 = Cos(y*Rad2Deg) ;Cos(y)
   If (cos2>0.0)
    scale1 = scale*Width/Height/cos2/Pi;
    Depth = 3*(Int(log_2(scale1*Height)))+3;
    For i = 0 To Width-1
     theta1 = longi-0.5*Pi+Pi*(2.0*i-Width)/Width/scale;
     col(i,j) = planet0(Cos(theta1*Rad2Deg)*cos2,Sin(y*Rad2Deg),-Sin(theta1*Rad2Deg)*cos2);
     ;col(i,j) = planet0(Cos(theta1)*cos2,Sin(y),-Sin(theta1)*cos2);
     If (doshade) shades(i,j) = shade;
    Next
   EndIf
  EndIf
 Next

 If (hgrid <> 0.0) ;draw horisontal gridlines
  theta0 = 0.0
  While theta0>-90.0
  theta0=theta0-hgrid
  theta1 = theta0
   While theta1<90.0
    theta1=theta1+hgrid
    y = DEG2RAD*theta1;
    j = Height/2+Int(0.5*y*Width*scale/Pi)+k;
    If (j>=0 And j<Height)
     For i = 0 To Width-1: col(i,j) = BLACK: Next;
    EndIf
   Wend
  Wend
 EndIf

 If (vgrid <> 0.0) ;draw vertical gridlines
  theta0 = 0.0
  While theta0>-360.0
   theta0=theta0-vgrid
   theta1 = theta0
   While theta1<360.0
    theta1=theta1+vgrid
    i = Int(0.5*Width*(1.0+scale*(DEG2RAD*theta1-longi)/Pi));
    If (i>=0 And i<Width)
     For j = max(0,Height/2-Int(0.25*Pi*Width*scale/Pi)+k) To min(Height,Height/2+Int(0.25*Pi*Width*scale/Pi)+k)-1
      col(i,j) = BLACK;
     Next
    EndIf
   Wend
  Wend
 EndIf

End Function

Function planet0(x#,y#,z#)
 ;x#,y#,z#

 Local alt#;
 Local colour;

 alt = planet1(x,y,z);

 If (altColors)

  Local snow# = 0.125;
  Local tree# = snow * 0.5;
  Local bare# = (tree + snow) / 2.0;
    
  If (latic)
   snow = snow - (0.13 * (y*y*y*y*y*y));
   bare = bare - (0.12 * (y*y*y*y*y*y));
   tree = tree - (0.11 * (y*y*y*y*y*y));
  EndIf

  If (alt > 0) ;Land
   If (alt > snow) ;Snow: White
    colour = WHITE;
   ElseIf (alt > bare) ;Snow: Grey - White
    colour = GREY0+Int((1+LAND1-GREY0) * (alt-bare)/(snow-bare));
    If (colour > LAND1) colour = LAND1;
   ElseIf (alt > tree) ;Bare: Brown - Grey
    colour = GREEN1+Int((1+GREY0-GREEN1) * (alt-tree)/(bare-tree));
    If (colour > GREY0) colour = GREY0;
   Else ;Green: Green - Brown
    colour = LAND0+Int((1+GREEN1-LAND0) * (alt)/(tree));
    If (colour > GREEN1) colour = GREEN1;
   EndIf
  Else ;Sea
   alt = alt/2;
   If (alt > snow) ;Snow: White
    colour = WHITE;
   ElseIf (alt > bare)
    colour = GREY0+Int((1+LAND1-GREY0) * (alt-bare)/(snow-bare));
    If (colour > LAND1) colour = LAND1;
   Else
    colour = BLUE1+Int((BLUE1-BLUE0+1)*(25*alt));
    If (colour<BLUE0) colour = BLUE0;
   EndIf
  EndIf

 Else ;calculate colour

  If (alt <=0.0) ;if below sea level then
   If (latic And y*y+alt >= 0.98)
    colour = LAND4 ;white if close to poles
   Else
    colour = BLUE1+Int((BLUE1-BLUE0+1)*(10*alt)) ;blue scale otherwise
    If (colour<BLUE0) colour = BLUE0;
   EndIf
  Else
   If (latic) alt = alt + 0.10204*y*y ;altitude adjusted with latitude
   If (alt >= 0.1) ;if high then
    colour = LAND4;
   Else ;else green to brown scale
    colour = LAND0+Int((LAND2-LAND0+1)*(10*alt));
    If (colour>LAND2) colour = LAND2;
   EndIf
  EndIf

 EndIf

 Return (colour);

End Function

Function planet#(a#,b#,c#,d#,as#,bs#,cs#,ds#,ax#,ay#,az#,bx#,by#,bz#,cx#,cy#,cz#,dx#,dy#,dz#,x#,y#,z#,level)
 ;a#,b#,c#,d# = altitudes of the 4 vertices
 ;as#,bs#,cs#,ds# = seeds of the 4 vertices
 ;ax#,ay#,az#,bx#,by#,bz#,cx#,cy#,cz#,dx#,dy#,dz# = vertex coordinates
 ;x#,y#,z# = goal point
 ;level = levels to go

 Local ssa#,ssb#,ssc#,ssd#,ssas#,ssbs#,sscs#,ssds#,ssax#,ssay#,ssaz#
 Local ssbx#,ssby#,ssbz#,sscx#,sscy#,sscz#,ssdx#,ssdy#,ssdz#;

 Local abx#,aby#,abz#,acx#,acy#,acz#,adx#,ady#,adz#;
 Local bcx#,bcy#,bcz#,bdx#,bdy#,bdz#,cdx#,cdy#,cdz#;
 Local lab#,lac#,lad#,lbc#,lbd#,lcd#;
 Local ex#,ey#,ez#,e#,es#,es1#,es2#,es3#;
 Local eax#,eay#,eaz#,epx#,epy#,epz#;
 Local ecx#,ecy#,ecz#,edx#,edy#,edz#;
 Local x1#,y1#,z1#,x2#,y2#,z2#,l1#,tmp#;

 If (level>0)

  If (level=11)
   ssa=a: ssb=b: ssc=c: ssd=d: ssas=as: ssbs=bs: sscs=cs: ssds=ds;
   ssax=ax: ssay=ay: ssaz=az: ssbx=bx: ssby=by: ssbz=bz;
   sscx=cx: sscy=cy: sscz=cz: ssdx=dx: ssdy=dy: ssdz=dz;
  EndIf

  abx = ax-bx: aby = ay-by: abz = az-bz;
  acx = ax-cx: acy = ay-cy: acz = az-cz;
  lab = abx*abx+aby*aby+abz*abz;
  lac = acx*acx+acy*acy+acz*acz;

  If (lab<lac)
   Return (planet(a,c,b,d,as,cs,bs,ds,ax,ay,az,cx,cy,cz,bx,by,bz,dx,dy,dz,x,y,z,level));
  Else
   adx = ax-dx: ady = ay-dy: adz = az-dz;
   lad = adx*adx+ady*ady+adz*adz;
   If (lab<lad)
    Return (planet(a,d,b,c,as,ds,bs,cs,ax,ay,az,dx,dy,dz,bx,by,bz,cx,cy,cz,x,y,z,level));
   Else
    bcx = bx-cx: bcy = by-cy: bcz = bz-cz;
    lbc = bcx*bcx+bcy*bcy+bcz*bcz;
    If (lab<lbc)
     Return (planet(b,c,a,d,bs,cs,as,ds,bx,by,bz,cx,cy,cz,ax,ay,az,dx,dy,dz,x,y,z,level));
    Else
     bdx = bx-dx: bdy = by-dy: bdz = bz-dz;
     lbd = bdx*bdx+bdy*bdy+bdz*bdz;
     If (lab<lbd)
      Return (planet(b,d,a,c,bs,ds,as,cs,bx,by,bz,dx,dy,dz,ax,ay,az,cx,cy,cz,x,y,z,level));
     Else
      cdx = cx-dx: cdy = cy-dy: cdz = cz-dz;
      lcd = cdx*cdx+cdy*cdy+cdz*cdz;
      If (lab<lcd)
       Return (planet(c,d,a,b,cs,ds,as,bs,cx,cy,cz,dx,dy,dz,ax,ay,az,bx,by,bz,x,y,z,level));
      Else
       es = rand2(as,bs);
       es1 = rand2(es,es);
       es2 = 0.5+0.1*rand2(es1,es1);
       es3 = 1.0-es2;
       If (ax=bx) ;very unlikely to ever happen
        ex = 0.5*ax+0.5*bx: ey = 0.5*ay+0.5*by: ez = 0.5*az+0.5*bz;
       ElseIf (ax<bx)
        ex = es2*ax+es3*bx: ey = es2*ay+es3*by: ez = es2*az+es3*bz;
       Else
        ex = es3*ax+es2*bx: ey = es3*ay+es2*by: ez = es3*az+es2*bz;
       EndIf
       If (lab>1.0) lab = lab^0.75 ;pow(lab,0.75);
       e = 0.5*(a+b)+es*dd1*Abs(a-b)+es1*dd2*(lab^POW) ;pow(lab,POW);
       eax = ax-ex: eay = ay-ey: eaz = az-ez;
       epx =  x-ex: epy =  y-ey: epz =  z-ez;
       ecx = cx-ex: ecy = cy-ey: ecz = cz-ez;
       edx = dx-ex: edy = dy-ey: edz = dz-ez;
       If ((eax*ecy*edz+eay*ecz*edx+eaz*ecx*edy-eaz*ecy*edx-eay*ecx*edz-eax*ecz*edy)*(epx*ecy*edz+epy*ecz*edx+epz*ecx*edy-epz*ecy*edx-epy*ecx*edz-epx*ecz*edy)>0.0)
        Return (planet(c,d,a,e,cs,ds,as,es,cx,cy,cz,dx,dy,dz,ax,ay,az,ex,ey,ez,x,y,z,level-1));
       Else
        Return (planet(c,d,b,e,cs,ds,bs,es,cx,cy,cz,dx,dy,dz,bx,by,bz,ex,ey,ez,x,y,z,level-1));
       EndIf
      EndIf
     EndIf
    EndIf
   EndIf 
  EndIf

 Else

  If (doshade)
   x1 = 0.25*(ax+bx+cx+dx);
   x1 = a*(x1-ax)+b*(x1-bx)+c*(x1-cx)+d*(x1-dx);
   y1 = 0.25*(ay+by+cy+dy);
   y1 = a*(y1-ay)+b*(y1-by)+c*(y1-cy)+d*(y1-dy);
   z1 = 0.25*(az+bz+cz+dz);
   z1 = a*(z1-az)+b*(z1-bz)+c*(z1-cz)+d*(z1-dz);
   l1 = Sqr(x1*x1+y1*y1+z1*z1);
   If (l1=0.0) l1 = 1.0;
   tmp = Sqr(1.0-y*y);
   If (tmp<0.0001) tmp = 0.0001;
   x2 = x*x1+y*y1+z*z1;
   y2 = -x*y/tmp*x1+tmp*y1-z*y/tmp*z1;
   z2 = -z/tmp*x1+x/tmp*z1;
   shade = Int((-Sin(Pi*shade_angle/180.0)*y2-Cos(Pi*shade_angle/180.0)*z2)/l1*48.0+128.0);
   If (shade<10) shade = 10;
   If (shade>255) shade = 255;
  EndIf
  Return ((a+b+c+d)/4);

 EndIf

End Function

Function planet1#(x#,y#,z#)
 ;x#,y#,z#

 Local abx#,aby#,abz#,acx#,acy#,acz#,adx#,ady#,adz#,apx#,apy#,apz#;
 Local bax#,bay#,baz#,bcx#,bcy#,bcz#,bdx#,bdy#,bdz#,bpx#,bpy#,bpz#;

 abx = ssbx-ssax: aby = ssby-ssay: abz = ssbz-ssaz;
 acx = sscx-ssax: acy = sscy-ssay: acz = sscz-ssaz;
 adx = ssdx-ssax: ady = ssdy-ssay: adz = ssdz-ssaz;
 apx = x-ssax: apy = y-ssay: apz = z-ssaz;

 If ((adx*aby*acz+ady*abz*acx+adz*abx*acy-adz*aby*acx-ady*abx*acz-adx*abz*acy)*(apx*aby*acz+apy*abz*acx+apz*abx*acy-apz*aby*acx-apy*abx*acz-apx*abz*acy)>0.0);
  ;p is on same side of abc as d
  If ((acx*aby*adz+acy*abz*adx+acz*abx*ady-acz*aby*adx-acy*abx*adz-acx*abz*ady)*(apx*aby*adz+apy*abz*adx+apz*abx*ady-apz*aby*adx-apy*abx*adz-apx*abz*ady)>0.0);
   ;p is on same side of abd as c
   If ((abx*ady*acz+aby*adz*acx+abz*adx*acy-abz*ady*acx-aby*adx*acz-abx*adz*acy)*(apx*ady*acz+apy*adz*acx+apz*adx*acy-apz*ady*acx-apy*adx*acz-apx*adz*acy)>0.0);
    ;p is on same side of acd as b
    bax = -abx: bay = -aby: baz = -abz;
    bcx = sscx-ssbx: bcy = sscy-ssby: bcz = sscz-ssbz;
    bdx = ssdx-ssbx: bdy = ssdy-ssby: bdz = ssdz-ssbz;
    bpx = x-ssbx: bpy = y-ssby: bpz = z-ssbz;
    If ((bax*bcy*bdz+bay*bcz*bdx+baz*bcx*bdy-baz*bcy*bdx-bay*bcx*bdz-bax*bcz*bdy)*(bpx*bcy*bdz+bpy*bcz*bdx+bpz*bcx*bdy-bpz*bcy*bdx-bpy*bcx*bdz-bpx*bcz*bdy)>0.0);
     ;p is on same side of bcd as a
     ;Hence, p is inside tetrahedron
     Return (planet(ssa,ssb,ssc,ssd,ssas,ssbs,sscs,ssds,ssax,ssay,ssaz,ssbx,ssby,ssbz,sscx,sscy,sscz,ssdx,ssdy,ssdz,x,y,z,11));
    EndIf
   EndIf
  EndIf
 EndIf

 ;otherwise
 Return (planet(M,M,M,M,r1,r2,r3,r4,0.0,0.0,3.01,0.0,Sqr(8.0)+0.01*r1*r1,-1.02+0.01*r2*r3,-Sqr(6.0)-0.01*r3*r3,-Sqr(2.0)-0.01*r4*r4,-1.02+0.01*r1*r2,Sqr(6.0)-0.01*r2*r2,-Sqr(2.0)-0.01*r3*r3,-1.02+0.01*r1*r3,x,y,z,Depth));
 ;M,M,M,M = initial altitude is M on all corners of tetrahedron
 ;r1,r2,r3,r4 = same seed set is used in every call
 ;coordinates of vertices
 ;x,y,z = coordinates of point we want colour of
 ;Depth = subdivision depth

End Function

Function rand2#(p#,q#)
 ;p#,q#
 ;random number generator taking two seeds
 ;rand2(p,q) = rand2(q,p) is important

 Local r#;
 r = (p+3.14159265)*(q+3.14159265);
 Return (2.0*(r-Int(r))-1.0);

End Function

Function log_2#(x#)
 Return (Log(x)/Log(2.0))
End Function

Subirenihil

2006

Top row of pixels is still black. Otherwise VERY nice. :)

puki

2006

Wow - I totally love this.

I'm claiming the code.

Mine, mine, mine.

Pineapple

2006

Great work!

Dabz

Code Archives Forum