#include "functions.inc"

#declare rds=false;

global_settings
{
  assumed_gamma 2.2
  #if (rds = true)
    #declare amb=0;
    radiosity { } // most basic radiosity, to be detailed in the future!
  #else
    #declare amb=0.05;
  #end
}


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


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

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

#declare latstart=41;
#declare latend=42;
#declare longstart=44;
#declare longend=45;
#declare latsize=latend-latstart;
#declare longsize=longend-longstart;
#declare month=7;

#declare res=36; // 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 NumTextures = NumVertices;
#declare V_vec_Arr = array[NumVertices];
#declare Face_Arr = array[NumFaces];
#declare N_vec_Arr = array[NumNormals];
#declare T_Arr = array[NumTextures];

#declare a = 0;

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

#declare P_Texture=
pigment
{
  image_map
  {
    "Texturen/n41e044_2004-07.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))*(1-flattening)*(radians(latstart+(a/ydim)*latsize)),
		(1-flattening)*sin(radians(latstart+(a/ydim)*latsize)), cos(radians(-longstart-(b/xdim)*longsize))*(1-flattening)*cos(radians(latstart+(a/ydim)*latsize))>;
    #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 F_Earthslice = finish {ambient ",str(amb, 1, 2)," diffuse 1 brilliance 0.25 }\n"))
#write (ES, concat("#declare Earth_Slice_",tilename,"=\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], ",",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")
#if (!(tilename="n80e014" | tilename="n80e016" | tilename="n80e017"))
  #declare i=0;
  #write (ES, "  texture_list")
  #write (ES, "  {\n")
  #write (ES, concat("    ",str(NumTextures, 1, 0),"\n"))
  #declare a=0;
  #declare c=0;
  #while (a<ydim)
    #declare b=0;
    #while (b<xdim)
      #declare C_Texture = eval_pigment(P_Texture, <(0.5+b)*(1/xdim), (0.5+a)*(1/ydim), 1>);
      #write (ES, "    texture\n")
      #write (ES, "    {\n")
      #write (ES, "      pigment\n")
      #write (ES, "      {\n")
      #write (ES, concat("	   color rgb <",vstr(3, C_Texture, ",", 1, 7),">\n"))
      #write (ES, "      }\n")
      #write (ES  "	 finish{ F_Earthslice }\n")
      #write (ES, "    }\n"
      #declare c=c+1;
      #declare b=b+1;
    #end
    #declare a=a+1;
  #end
  #write (ES, "  }\n")
#end


#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")




/* #while (i < NumTextures)
  #write (ES, concat("texture { T_Arr[",str(i, 1, 0)))
  #if (i=NumTextures-1)
    #write (ES, "\n")  
  #else
    #write (ES, ",\n")
  #end
  #declare i=i+1; 
#end */

#if (tilename="n80e014" | tilename="n80e016" | tilename="n80e017")
  #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, "  double_illuminate\n")
/* #write (ES, "  translate ApproxLook\n") */
#write (ES, "}\n")
#fclose ES