// Scene file for mapping an astronomical image defined
// in galactic coordinates onto the celestial sphere for
// an observer on Earth at a given time at location.
//
// Copyright (C) 2009 by Christian Froeschlin, www.chrfr.de
//
// Version 1.0
//
// ESO image http://www.gigagalaxyzoom.org/B.html created by
// Serge Brunier, brought to my attention by a post from "waggy".
// Reuses code from Ingo Janssen's SunPos macro. The online tool
// http://violet.pha.jhu.edu/cgi-bin/eqtogal_tool was used to
// obtain coordinates for vernal point and equatorial north
// pole in galactic coordinates.
//
// Issues: 1. Precession not taken into account
//         2. The celestial equator in constellations.jpg
//            did not align perfectly with the horizon after
//            converting to equatorial coordinates, differing
//            about 2 degrees. Not sure why.


global_settings 
{
  assumed_gamma 1.0
}


// Configure observer time & location

#declare LONG = 15;
#declare LAT  = 50;

#declare YEAR  = 2000;
#declare MONTH =   03;
#declare DAY   =   20;

#declare HOUR  =   06;
#declare MIN   =   00;


// Configure viewing direction

camera
{
  location <0,0,0>
  //look_at z  // North
  look_at x  // Due East
  //look_at y  // Up
  angle 90
}


// Configure galaxy if godlike powers available, otherwise don't modify

// Vernal point and equatorial north pole in galactic coordinates.
// The values actually depend on precession but should be sufficiently
// precise for visualization purposes between 1500 and 2500 AD.
#declare L_VERNAL =  96.34;
#declare B_VERNAL = -60.19;
#declare L_POLE   = 302.93;
#declare B_POLE   = -27.13;



// The code to calculate sidereal time from local time for a given location
// was already available in the SunPos macro. Here we basically define a direct
// copy of that macro with a different return value.
#macro HourAngle(Year, Month, Day, Hour, Minute, Lstm, LAT, LONG)
   
   // Calculate universal time (UT)
   #local T= Hour+(Minute/60);
   #local UT= T-Lstm/15;
   #if (0>UT)
       #local Day= Day-1;
       #local UT= 24+UT;
   #end
   #if (UT>24)
       #local Day= Day+1;
       #local UT= UT-24;
   #end
   
   // Amount of days to, or from, the year 2000
   #local d= 367*Year-int((7*int((Year+int((Month+9))/12)))/4)+int((275*Month)/9)+Day-730530+UT/24;
   
   // Longitude of perihelion (w), eccentricity (e)
   #local w= 282.9404+4.70935E-5*d;
   #local e= 0.016709-1.151E-9*d;
   
   // Mean anomaly (M), sun's mean longitude (L)
   #local M= 356.0470+0.9856002585*d;
   #if (0<M<360)
       #local M= M-floor(M/360)*360;
   #end
   #local L= w+M;
   #if (0<L<360)
       #local L= L-floor(L/360)*360;
   #end
   
   // Obliquity of the ecliptic, eccentric anomaly (E)
   #local oblecl= 23.4393-3.563E-7*d;
   #local E= M+(180/pi)*e*sin(radians(M))*(1+e*cos(radians(M)));
   
   // Sun's rectangular coordinates in the plane of ecliptic (A,B)
   #local A= cos(radians(E))-e;
   #local B= sin(radians(E))*sqrt(1-e*e);
   
   // Distance (r), true anomaly (V), longitude of the sun (lon)
   #local r= sqrt(A*A+B*B);
   #local V= degrees(atan2(radians(B),radians(A)));
   #local lon= V+w;
   #if (0<lon<360)
       #local lon= lon-floor(lon/360)*360;
   #end
   
   // Calculate declination (Decl) and right ascension (RA)
   #local Decl= degrees(asin(sin(radians(oblecl))*sin(radians(lon))));
   #local RA= degrees(atan2(sin(radians(lon))*cos(radians(oblecl)),cos(radians(lon))))/15;
   
   // Greenwich meridian siderial time at 00:00 (GMST0), siderial time (SIDTIME), hour angle (HA)
   #local GMST0= L/15+12;
   #local SIDTIME= GMST0+UT+LONG/15;
   #local HA = SIDTIME*15;
   HA
#end



// Actual star map
sphere 
{
  0, 10000
  pigment 
  {
    image_map
    {
      png "starmap.png" 
      //png "lowres.png"  
      //jpeg "constellations.jpg"  
      map_type    1     
      interpolate 4 // 0=none, 1=linear, 2=bilinear, normalized
      once
    }
    // Normalize spherical projection of image to right-handed
    // galactic cooordinate system with galactic center at +z.
    // Galactic plane remains in x-z-plane
    scale -1
    rotate -90*y 
  }
  finish {ambient 1}
  
      
  // *****************************************************
  // Transform to equatorial coordinate system with spring
  // equinox at +z and equatorial plane in x-z-plane. 
  // *****************************************************


  // Original location of vernal point (the coordinates could
  // be transformed directly but creating this point and using
  // vrotate avoids having to do the spherical trigonometry).
  // Also, the point can be used for debug visualization.
  // Note :The vernal point is located between the constellations
  // Pisces and Aquarius where the celestial pole and the
  // ecliptic intersects (displayed in constellations.jpg)
  #declare P_VERNAL = 100*z;
  #declare P_VERNAL = vrotate (P_VERNAL,B_VERNAL*x); 
  #declare P_VERNAL = vrotate (P_VERNAL,L_VERNAL*y);

  // First get the celestial north pole to +y
  rotate -L_POLE*y
  rotate (90-B_POLE)*x
    
  // This transformation also affects the vernal point:
  #declare P_VERNAL = vrotate (P_VERNAL,-L_POLE*y); 
  #declare P_VERNAL = vrotate (P_VERNAL,(90-B_POLE)*x);

  // Sanity checks:  P_VERNAL should now lie in the x-z plane
  //                 The celestial equator should now be on the "horizon"
  //                 Polaris should be at +y           

  // Now correct ascension such that the vernal point ends up at +z

  #declare A_VERNAL = 90-degrees(atan2(P_VERNAL.z,P_VERNAL.x)); 
  rotate -A_VERNAL*y

  // Sanity check: The vernal point should now be at +z.


  // *****************************************************
  // Transform to horzizon coordinate system with north
  // at +z and horizon in the x-z plane 
  // *****************************************************


  // By definition, the hour angle of the vernal point is the
  // sidereal time. 
  #declare A_LOCAL = HourAngle(YEAR,MONTH,DAY,HOUR,MIN,LONG,LAT,LONG);
  rotate -A_LOCAL*y
 
  // Polaris is still at +y, it's final position is determined by latitude
  rotate (90-LAT)*x

  // Sanity check: Around March 20th, the vernal point should be visible
  //               due east at 6am local time from any location excepting
  //               the poles.
}


// Ground plane
plane {y,-0.1 pigment {color rgbt <0,1,0,0.3>} finish {ambient 0.3}}