// Spherical mesh section from SRTM data
// Thanks to Jerome M. Berger for improving the algorithm and to Andrel Linnenbank to correct matrix size!



#declare l=1201;

#declare mtrix=array[l][l]

#declare oHgh=1.7; // Observer height above ground

#declare cLat=34.3282;    // Camera location latitude
#declare cLong=-68.5;  // Camera location longitude
#declare lLat=34.5782;    // Camera look_at latitude
#declare lLong=-68.0765;  // Camera look_at longitude

// 1 POV unit = 1 metre

#declare re=6378140; // Earth, equatorial radius
#declare rp=6356755; // Earth, polar radius

#declare obl=rp/re; // Earth oblateness factor

#declare sLat=34;
#declare eLat=35;
#declare sLong=-68;
#declare eLong=-69;

#declare os=clock-int(clock); // flag for operating system: 0 = Linux, 0.5 = Windows

#declare F_Standard=
finish { ambient 0.1 diffuse 1 brilliance 0.5 }

#warning "Started reading ASCII matrix"
    #switch(os)
      #case (0)
        #fopen matr "/windows/e/homepage/n34e068-b.txt" read
      #break
      #case (0.5)
        #fopen matr "e:\homepage\n34e068-b.txt" read
      #break
      #default
        #fopen matr "/windows/e/homepage/n34e068-b.txt" read
      #break
    #end
    #declare a = 0;
    #while (a < l)
      #declare b = 0;
      #while (b < l)
        #read (matr, data)
        #declare mtrix[a][b]=data;
        #declare b = b+1;
      #end
      #warning concat("Reading line ", str(a, 4, 0))
      #declare a = a+1;
    #end
    #fclose matr
#warning "Reading of ASCII matrix completed!"

#declare EarthSlice=
mesh2 {
  vertex_vectors 
  {
    l*l,
    #declare a=0;
    #while (a<l)
      #declare b=0;
      #while (b<l)
	<(re+mtrix[a][b]) * sin ((sLong-b/l) * (pi/180)) *
        cos ((sLat + a/l) * (pi/180)),
        (rp+mtrix[a][b])*sin((sLat+a/l)*(pi/180)),
        (re+mtrix[a][b]) * cos ((sLong-b/l) * (pi/180)) *
        cos ((sLat+a/l)*(pi/180))>
        #if (b<l-1)
          ,
        #end
        #declare b=b+1;
      #end
      #warning concat("Assigning vectors for line ", str(a, 4, 0))
      #declare a=a+1;
    #end
  }
  face_indices {
    (l-1) * (l-1) * 2,
    #declare a = 0;
    #while (a < l-1)
      #declare b = 0;
      #while (b < l-1)
        <l*a + b, l*a + b+1, l*(a+1) + b>,
        <l*a + b+1, l*(a+1) + b+1, l*(a+1) + b>
        #if(b<(l-2))
          ,
        #end
        #declare b = b+1;
      #end
      #warning concat("Parsing lines ", str(a, 4, 0), " and ", str(a+1, 4, 0),".")
      #declare a = a+1;
    #end
  }
  texture {
    pigment { color rgb <0.9, 0.81, 0.4> }
    finish { F_Standard }
  }
} 


object { EarthSlice }  

light_source
{
  <-5000000, 30000000, 10000000>
  color rgb 1
}

#declare geoVect=vnormalize(<sin(cLong*(pi/180))*cos(cLat*(pi/180)), sin(cLat*(pi/180)), cos(cLong*(pi/180))*cos(cLat*(pi/180))>);

#declare cNormVect=10000000*geoVect;
#declare lNormVect=10000000*vnormalize(<sin(lLong*(pi/180))*cos(lLat*(pi/180)), sin(lLat*(pi/180)), cos(lLong*(pi/180))*cos(lLat*(pi/180))>);

#declare lo=trace(EarthSlice, cNormVect, -geoVect)+geoVect*1.7;
#declare la=trace(EarthSlice, lNormVect, -geoVect);

camera
{
  angle 40
  sky geoVect
  location lo
  look_at la
  
  // lNormVect/10000000*vlength(trace(EarthSlice, lNormVect, -geoVect))
}

sky_sphere
{
  pigment { color rgb <0.3, 0.6, 0.8> }
}