// Spherical mesh section from SRTM data
// Thanks to Jerome M. Berger for improving the algorithm!

#declare mtrix=array[1200][1200]



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

#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 < 1200)
      #declare b = 0;
      #while (b < 1200)
        #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!"

mesh2 {
  vertex_vectors 
  {
    1200*1200,
    #declare a=0;
    #while (a<1200)
      #declare b=0;
      #while (b<1200)
	<(re+mtrix[a][b]) * sin ((sLong-b/1200) * (pi/180)) *
        cos ((sLat + a/1200) * (pi/180)),
        (rp+mtrix[a][b])*sin((sLat+a/1200)*(pi/180)),
        (re+mtrix[a][b]) * cos ((sLong-b/1200) * (pi/180)) *
        cos ((sLat+a/1200)*(pi/180))>
        #if (b<1199)
          ,
        #end
        #declare b=b+1;
      #end
      #warning concat("Assigning vectors for line ", str(a, 4, 0))
      #declare a=a+1;
    #end
  }
  face_indices {
    1199 * 1199 * 2,
    #declare a = 0;
    #while (a < 1199)
      #declare b = 0;
      #while (b < 1199)
        <1200*a + b, 1200*a + b+1, 1200*(a+1) + b>,
        <1200*a + b+1, 1200*(a+1) + b+1, 1200*(a+1) + b>
        #if(b<1198)
          ,
        #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 <1, 0.5, 0> }
    finish { F_Standard }
  }
} 


sphere
{
  0, re+2500
  scale <1, rp/re, 1>
  texture
  {
    pigment { color rgb <0, 0, 1> }
    finish { F_Standard }
  }
} 
  

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

camera
{
  location 6800000*<sin(-68.5*(pi/180))*cos(34.5*(pi/180)), sin(34.5*(pi/180)), cos(-68.5*(pi/180))*cos(34.5*(pi/180))>
  look_at 0
  angle 40
}

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