/** This generates a model of the moon for 3-D printing.

    It requires data files from NASA:
    https://svs.gsfc.nasa.gov/4720

    Specifically the 1440x720 16-bit unsigned int TIFF displacement file
    renamed here to "moonhalfmeter.tif".  You can also use larger files but
    you'll have to scale them down by setting the "rescale" parameter below.

    Use the solver
https://frinklang.org/fsp/solve2.fsp?equations=%2F%2F+Globe+solver+for+3-D+printing%0D%0Ad+%3D+2+r%0D%0Ac+%3D+pi+d%0D%0Apixelsize+%3D+width+%2F+c%0D%0Apixelsize+%3D+res%0D%0A%0D%0Ascale+%3D+r+%2F+moonradius%0D%0A&solveFor=&ev=on&sel_c=S&val_c=299792458+m+s%5E-1&sel_d=S&val_d=86400+s&sel_moonradius=S&val_moonradius=1.738000e%2B6+m&sel_pi=L&val_pi=3.141592653589793238&sel_pixelsize=S&val_pixelsize=&sel_r=S&val_r=&sel_res=L&val_res=254+pixels%2Fin&sel_scale=S&val_scale=&sel_width=L&val_width=1440+pixels&resultAs=&showorig=on

*/

res = 254 / in   // Resolution of the model in voxels/inch (0.1 mm)

resample = 1  // This is the factor to scale down a very large image.  Use 1
              // to use full resolution.  Use 2 to sample every other pixel.

img = new image["file:moonhalfmeter.tif"]
[width, height] = img.getSize[]

hmax = 0
hmin = 255

for x = 0 to width-1
   for y = 0 to height-1
   {
      h = img.getPixelGrayInt[x,y]
      if h < hmin
         hmin = h
      if h > hmax
         hmax = h
   }

println["hmin is $hmin, hmax is $hmax"]

// exaggeration factor of vertical scale
exaggeration = 10

// Wall thickness of model.
thickness = 1.19 mm   
   
/* According to the page, the unsigned int tiff files are in units of half-
   meters, relative to a radius of 1727400 m. */
filebaseradius = 1727400 m

moonMeanRadius = 1737.4 km   // LRO reference sphere

maxMoonRadius = filebaseradius + hmax * 1/2 m * (2^16 / 256) * exaggeration
minMoonRadius = filebaseradius + hmin * 1/2 m * (2^16 / 256) * exaggeration

scaledWidth = floor[width/resample]
scaledHeight = floor[height/resample]

r = scaledWidth / (2 pi res)
println["r = $r"]

scale = r / maxMoonRadius
println["scale = $scale"]

println["r res = " + (r res)]
v = callJava["frink.graphics.VoxelArray", "construct", [-r res, r res, -r res, r res, -r res, r res, false]]

baseRadius = ceil[minMoonRadius scale res]
println["minMoonRadius scale res = " + baseRadius ]
baseSphere = callJava["frink.graphics.VoxelArray", "makeSphere", [baseRadius]]

v.add[baseSphere]   

// Make a 1-pixel tool   
tool = newJava["frink.graphics.VoxelArray", [1,1,1,true]]
   
for x = 0 to (width-1) step resample
{
   sx = x / resample
   long = (x circle / width)
   clong = cos[long]
   slong = sin[long]

   for y = 0 to (height-1) step resample
   {
      sy = y / resample
      lat = ((height - y) / height) 180 deg - 90 deg
      clat = cos[lat]
      slat = sin[lat]
      
      h = img.getPixelGrayInt[x,y]
      rh = (filebaseradius + h * 1/2 m * (2^16 / 256) * exaggeration) scale
      highX = rh clat clong
      highY = rh clat slong
      highZ = rh slat

      rl = baseRadius
      lowX = rl clat clong
      lowY = rl clat slong
      lowZ = rl slat

//      println["rl is $rl, rh is " + (rh res)]

      v.addAlongLine[tool, round[lowX], round[lowY], round[lowZ],
                     round[highX res], round[highY res], round[highZ res],
                     0, 0, 0]
   }
}

v.projectX[undef].show["X"]
v.projectY[undef].show["Y"]
v.projectZ[undef].show["Z"]

filename = "moon3DSolid.obj"
print["Writing $filename..."]
w = new Writer[filename]
w.println[v.toObjFormat["moon3dsolid", 1/(res mm)]]
w.close[]
println["done."]