#include "functions.inc"

global_settings
{
  assumed_gamma 2.2
}

#declare path="/home/Raytracing/Povscn/POVEarth/Rohdaten/"; // path for ASCII file containg elevation data
#declare path_textures="/home/Raytracing/Povscn/POVEarth/Texturen/";    

#declare xdim=3601; // full resolution: 3601;
#declare ydim=3601; // full resolution: 3601;

#declare rd=6378140-422;
#declare rng=9270;

#declare latstart=80;
#declare latend=81;
#declare longstart=14;
#declare longend=15;
#declare latsize=latend-latstart;
#declare longsize=longend-longstart;
#declare month=7;

#declare res=3601; // full resolution: 3601;

#if (latstart-1 < 0)
  #declare lat="s";
#else
  #declare lat="n";
#end

#if (longstart < 0)
  #declare long="w";
#else
  #declare long="e";
#end

#declare tilename=concat(lat, str(latstart, -2, 0), long, str(longstart, -3, 0));

// tiles are named by the coordinates of their southwestern corner!

#declare altdataname = concat(path, tilename, ".txt");

#declare NumVertices = xdim*ydim;
#declare NumFaces = (xdim-1)*(ydim-1)*2;
#declare NumNormals = NumVertices;
#declare V_vec_Arr = array[NumVertices];
#declare Face_Arr = array[NumFaces];
#declare N_vec_Arr = array[NumNormals];

#declare a = 0;

#declare P_Heightfield=
pigment
{
  image_map
  {
    "Heightfields/n80e014_cleaned.png"
    interpolate 2
  }
}

#declare a=0;
#declare c=0;
#while (a < ydim)
  #declare b=0;
  #while (b < xdim)
    #declare hcolor = eval_pigment(P_Heightfield, <(0.5+b)*(1/xdim), (0.5+a)*(1/ydim), 1>);
    #declare hval = hcolor.red*rng+hcolor.green*(rng/255);
    #declare V_vec_Arr[c] = (rd+hval)*<sin(radians(-longstart-(b/xdim)*longsize))*cos(radians(latstart+(a/ydim)*latsize)),
		sin(radians(latstart+(a/ydim)*latsize)), cos(radians(-longstart-(b/xdim)*longsize))*cos(radians(latstart+(a/ydim)*latsize))>/* -ApproxLook */;
    #declare c=c+1;		  
    #declare b=b+1;
  #end
  #declare a=a+1;
#end

#declare a = 0;
#declare c = 0;
#while (a < ydim-1)
  #declare b = 0;
  #while (b < xdim-1)
    #declare Face_Arr[c] = <a*res+b, a*res+b+1, (a+1)*res+b>;
    #declare Face_Arr[c+1] = <a*res+b+1, (a+1)*res+b+1, (a+1)*res+b>;
    #declare c = c+2;
    #declare b = b+1;
  #end
  #declare a = a+1;
#end

#declare a=0; // vcross(C-B,A-B)
#declare c=0;
#while (a < ydim)
  #declare b=0;
  #while (b < xdim)
    #if (a=0 & b=0) // lower left corner: 1 adjacent face
      #declare nno = vcross(V_vec_Arr[xdim]-V_vec_Arr[1], V_vec_Arr[0]-V_vec_Arr[1]);
    #else
      #if (a=0 & b>0 & b<xdim-1) // row 0 except lower left and lower right corner: 3 adjacent faces per vertex
	#declare nno = vcross(V_vec_Arr[(a+1)*xdim+(b-1)]-V_vec_Arr[b], V_vec_Arr[b-1]-V_vec_Arr[b]);
	#declare nno = nno + vcross(V_vec_Arr[(a+1)*ydim+(b-1)]-V_vec_Arr[(a+1)*xdim+b], V_vec_Arr[b]-V_vec_Arr[(a+1)*xdim+b]);
	#declare nno = nno + vcross(V_vec_Arr[(a+1)*ydim+b]-V_vec_Arr[(a+1)*xdim+b+1], V_vec_Arr[b]-V_vec_Arr[(a+1)*xdim+b+1]);
	#declare nno = nno/3;
      #else
	#if (a=0 & b=xdim-1) // lower right corner: 2 adjacent faces
	  #declare nno = vcross(V_vec_Arr[xdim*2-2]-V_vec_Arr[xdim-1], V_vec_Arr[xdim-2]-V_vec_Arr[xdim-1]);
	  #declare nno = nno + vcross(V_vec_Arr[xdim*2-2]-V_vec_Arr[xdim*2-1], V_vec_Arr[xdim-1]-V_vec_Arr[xdim*2-1]);
	  #declare nno = nno/2;
	#else
	  #if (a>0 & a<ydim-1 & b=0) // leftmost column except lower left and upper left corner: 3 adjacent faces per vertex
	    #declare nno = vcross(V_vec_Arr[a*ydim]-V_vec_Arr[(a-1)*ydim+1], V_vec_Arr[(a-1)*ydim]-V_vec_Arr[(a-1)*ydim+1]);
	    #declare nno = nno + vcross(V_vec_Arr[a*ydim+1]-V_vec_Arr[(a-1)*ydim+1], V_vec_Arr[a*ydim]-V_vec_Arr[(a-1)*ydim+1]);
	    #declare nno = nno + vcross(V_vec_Arr[(a+1)*ydim]-V_vec_Arr[a*ydim+1], V_vec_Arr[a*ydim]-V_vec_Arr[a*ydim+1]);
	    #declare nno = nno/3;
	  #else
	    #if (a>0 & a<ydim-1 & b>0 & b<xdim-1) // all rows and columns not touching any edge: 6 adjacent faces per vertex
	      #declare nno = vcross(V_vec_Arr[(a-1)*xdim+b]-V_vec_Arr[a*xdim+b-1], V_vec_Arr[a*xdim+b]-V_vec_Arr[a*xdim+b-1]);
	      #declare nno = nno + vcross(V_vec_Arr[(a-1)*xdim+b+1]-V_vec_Arr[(a-1)*xdim+b], V_vec_Arr[a*xdim+b]-V_vec_Arr[(a-1)*xdim+b]);
	      #declare nno = nno + vcross(V_vec_Arr[a*xdim+b+1]-V_vec_Arr[(a-1)*xdim+b+1], V_vec_Arr[a*xdim+b]-V_vec_Arr[(a-1)*xdim+b+1]);
	      #declare nno = nno + vcross(V_vec_Arr[(a+1)*xdim+b]-V_vec_Arr[a*xdim+b+1], V_vec_Arr[a*xdim+b]-V_vec_Arr[a*xdim+b+1]);
	      #declare nno = nno + vcross(V_vec_Arr[(a+1)*xdim+b-1]-V_vec_Arr[(a+1)*xdim+b], V_vec_Arr[a*xdim+b]-V_vec_Arr[(a+1)*xdim+b]); 
	      #declare nno = nno + vcross(V_vec_Arr[a*xdim+b-1]-V_vec_Arr[(a+1)*xdim+b-1], V_vec_Arr[a*ydim+b]-V_vec_Arr[(a+1)*xdim+b-1]);
	      #declare nno = nno/6;
	    #else
	      #if (a>0 & a<ydim-1 & b=ydim-1) // rightmost column except lower right and upper right corner: 3 adjacent faces per vertex
		#declare nno = vcross(V_vec_Arr[(a-1)*xdim+(xdim-1)]-V_vec_Arr[a*xdim+(xdim-2)], V_vec_Arr[a*xdim+(xdim-1)]-V_vec_Arr[a*xdim+(xdim-2)]);
		#declare nno = nno + vcross(V_vec_Arr[a*xdim+(xdim-2)]-V_vec_Arr[(a+1)*xdim+(xdim-2)], V_vec_Arr[a*xdim+(xdim-1)]-V_vec_Arr[(a+1)*xdim+(xdim-2)]);
		#declare nno = nno + vcross(V_vec_Arr[(a+1)*xdim+(xdim-2)]-V_vec_Arr[(a+1)*xdim+(xdim-1)], V_vec_Arr[a*xdim+(xdim-1)]-V_vec_Arr[(a+1)*xdim+(xdim-1)]);
		#declare nno = nno/3;
	      #else
		#if (a=ydim-1 & b=0) // upper left corner: 2 adjacent faces
		  #declare nno = vcross(V_vec_Arr[(a-1)*xdim+1]-V_vec_Arr[(a-1)*xdim], V_vec_Arr[a*xdim]-V_vec_Arr[(a-1)*xdim]);
		  #declare nno = nno + vcross(V_vec_Arr[a*xdim+1]-V_vec_Arr[(a-1)*xdim+1], V_vec_Arr[a*xdim]-V_vec_Arr[(a-1)*xdim+1]);
		  #declare nno = nno/2;
		#else
		  #if (a=ydim-1 & b>0 & b<xdim-1) // row 3600 except upper left and upper right corner: 3 adjacent faces per vertex
		    #declare nno = vcross(V_vec_Arr[(a-1)*xdim+b]-V_vec_Arr[a*xdim+b-1], V_vec_Arr[a*xdim+b]-V_vec_Arr[a*xdim+b-1]);
		    #declare nno = nno + vcross(V_vec_Arr[(a-1)*xdim+b+1]-V_vec_Arr[(a-1)*xdim+b], V_vec_Arr[a*xdim+b]-V_vec_Arr[(a-1)*xdim+b]);
		    #declare nno = nno + vcross(V_vec_Arr[a*xdim+b+1]-V_vec_Arr[(a-1)*xdim+b+1], V_vec_Arr[a*xdim+b]-V_vec_Arr[(a-1)*ydim+b+1]);
		    #declare nno = nno/3;
		  #else
		    #if (a=ydim-1 & b=xdim-1) // upper right corner, 1 adjacent face
		      #declare nno = vcross(V_vec_Arr[xdim*ydim-xdim]-V_vec_Arr[xdim*ydim-2], V_vec_Arr[xdim*ydim-1]-V_vec_Arr[xdim*ydim-2]);
		    #end
		  #end
		#end
	      #end
	    #end
	  #end
	#end
      #end
    #end
    #warning str(c, 1, 0)
    #declare N_vec_Arr[c] = vnormalize(nno);
    #declare b=b+1;
    #declare c=c+1;
  #end
  #declare a=a+1;
#end


#declare lNormVect=vnormalize(<sin(radians(-longstart-0.5))*cos(radians(latstart+0.5)), sin(radians(latstart+0.5)), 
cos(radians(-longstart-0.5))*cos(radians(latstart+0.5))>);
#declare ApproxLook = lNormVect * rd;

#fopen ES concat(tilename, ".inc") write
#write (ES, concat("#declare ApproxLook = <", vstr(3, ApproxLook,",",3,4), ">;\n")) 
#write (ES, concat("#declare Earth_Slice =\nmesh2\n{\n  vertex_vectors\n     {\n",str(xdim*ydim, 1, 0),"\n"))

#declare sf=1;

#declare i=0;
#while (i<NumVertices)
  #write (ES, concat("    <",vstr(3, V_vec_Arr[i]-ApproxLook, ",",3,4),">"))
  #if (i=NumVertices-1)
    #write (ES, "\n")
  #else
    #write (ES, ",\n")
  #end
  #declare i=i+1;
#end
#write (ES,"    }\n")
#write (ES,"	normal_vectors\n")
#write (ES,"	{\n")
#write (ES, concat("    ",str(xdim*ydim, 1, 0),"\n"))
#declare i=0;
#while (i<NumNormals)
  #write (ES, concat("    <",vstr(3, N_vec_Arr[i], ",",3,4),">"))
  #if (i=NumNormals-1)
    #write (ES, "\n")
  #else
    #write (ES, ",\n")
  #end
  #declare i=i+1;
#end
#write (ES,"    }\n")
#write (ES, concat("  face_indices\n{\n    ",str((xdim-1)*2*(ydim-1), 7, 0),"\n"))

#declare i=0;
#while (i<NumFaces)
  #write (ES, concat("    <",vstr(3, Face_Arr[i], ",",0, 0),">"))
  #if (i=NumFaces-1)
    #write (ES, "\n")  
  #else
    #write (ES, ",\n")
  #end
  #declare i=i+1;
#end
#write (ES, "  }\n")
#write (ES, "  double_illuminate\n")
#write (ES, "  translate ApproxLook\n")
#if (tilename="n80e014")
  #write (ES, "  texture\n")
  #write (ES, "  {\n")
  #write (ES, "    pigment { color rgb <0.66, 0.65, 0.65> }\n")
  #write (ES, "    finish { ambient 0.1 diffuse 1 brilliance 0.4 }\n")
  #write (ES, "  }\n")
#end
#write (ES, "}")

#fclose ES