contrib/gel/mrc/vpgl/bgeo/bgeo_sun

Go to the documentation of this file.
00001 #include "bgeo_sun_pos.h"
00002 #include <vnl/vnl_math.h>
00003 #include <vcl_cmath.h>
00004 
00005 // Declaration of some constants
00006 static double twopi  = 2.0*vnl_math::pi;
00007 static double rad    = vnl_math::pi/180.0;
00008 static double invrad = 180.0/vnl_math::pi;
00009 static double dEarthMeanRadius  =  6371.01;   // In km
00010 static double dAstronomicalUnit =  149597890; // In km
00011 
00012 void bgeo_sun_pos(int year, int month, int day,
00013                   int hours, int minutes, int seconds,
00014                   double longitude, double latitude,
00015                   double& sun_azimuth, double& sun_elevation)
00016 {
00017 // Main variables
00018   double dElapsedJulianDays;
00019   double dDecimalHours;
00020   double dEclipticLongitude;
00021   double dEclipticObliquity;
00022   double dRightAscension;
00023   double dDeclination;
00024 
00025   // Auxiliary variables
00026   double dY;
00027   double dX;
00028 
00029   // Calculate difference in days between the current Julian Day
00030   // and JD 2451545.0, which is noon 1 January 2000 Universal Time
00031   {
00032     double dJulianDate;
00033     long int liAux1;
00034     long int liAux2;
00035     // Calculate time of the day in UT decimal hours
00036     dDecimalHours = hours + (minutes + seconds / 60.0 ) / 60.0;
00037     // Calculate current Julian Day
00038     liAux1 =(month-14)/12;
00039     liAux2=(1461*(year + 4800 + liAux1))/4 + (367*(month
00040             - 2-12*liAux1))/12- (3*((year + 4900
00041           + liAux1)/100))/4+day-32075;
00042     dJulianDate=(double)(liAux2)-0.5+hours/24.0;
00043     // Calculate difference between current Julian Day and JD 2451545.0
00044     dElapsedJulianDays = dJulianDate-2451545.0;
00045   }
00046 
00047   // Calculate ecliptic coordinates (ecliptic longitude and obliquity of the
00048   // ecliptic in radians but without limiting the angle to be less than 2*Pi
00049   // (i.e., the result may be greater than 2*Pi)
00050   {
00051     double dMeanLongitude;
00052     double dMeanAnomaly;
00053     double dOmega;
00054     dOmega=2.1429-0.0010394594*dElapsedJulianDays;
00055     dMeanLongitude = 4.8950630+ 0.017202791698*dElapsedJulianDays; // Radians
00056     dMeanAnomaly = 6.2400600+ 0.0172019699*dElapsedJulianDays;
00057     dEclipticLongitude = dMeanLongitude + 0.03341607*vcl_sin( dMeanAnomaly )
00058                        + 0.00034894*vcl_sin( 2*dMeanAnomaly )-0.0001134
00059                        - 0.0000203*vcl_sin(dOmega);
00060     dEclipticObliquity = 0.4090928 - 6.2140e-9*dElapsedJulianDays
00061                        + 0.0000396*vcl_cos(dOmega);
00062   }
00063 
00064   // Calculate celestial coordinates ( right ascension and declination ) in radians
00065   // but without limiting the angle to be less than 2*Pi (i.e., the result may be
00066   // greater than 2*Pi)
00067   {
00068     double dSin_EclipticLongitude;
00069     dSin_EclipticLongitude= vcl_sin( dEclipticLongitude );
00070     dY = vcl_cos( dEclipticObliquity ) * dSin_EclipticLongitude;
00071     dX = vcl_cos( dEclipticLongitude );
00072     dRightAscension = vcl_atan2( dY,dX );
00073     if ( dRightAscension < 0.0 ) dRightAscension = dRightAscension + twopi;
00074     dDeclination = vcl_asin( vcl_sin( dEclipticObliquity )*dSin_EclipticLongitude );
00075   }
00076 
00077   // Calculate local coordinates ( azimuth and zenith angle ) in degrees
00078   {
00079     double dGreenwichMeanSiderealTime;
00080     double dLocalMeanSiderealTime;
00081     double dLatitudeInRadians;
00082     double dHourAngle;
00083     double dCos_Latitude;
00084     double dSin_Latitude;
00085     double dCos_HourAngle;
00086     double dParallax;
00087     dGreenwichMeanSiderealTime = 6.6974243242 +
00088                                  0.0657098283*dElapsedJulianDays + dDecimalHours;
00089     dLocalMeanSiderealTime = (dGreenwichMeanSiderealTime*15 + longitude)*rad;
00090     dHourAngle = dLocalMeanSiderealTime - dRightAscension;
00091     dLatitudeInRadians = latitude*rad;
00092     dCos_Latitude = vcl_cos(dLatitudeInRadians);
00093     dSin_Latitude = vcl_sin(dLatitudeInRadians);
00094     dCos_HourAngle= vcl_cos(dHourAngle);
00095     sun_elevation = vcl_acos( dCos_Latitude*dCos_HourAngle*vcl_cos(dDeclination)
00096                               + vcl_sin(dDeclination)*dSin_Latitude);
00097     dY = -vcl_sin(dHourAngle);
00098     dX = vcl_tan(dDeclination)*dCos_Latitude - dSin_Latitude*dCos_HourAngle;
00099     sun_azimuth = vcl_atan2( dY, dX );
00100     if ( sun_azimuth < 0.0 )
00101       sun_azimuth = sun_azimuth + twopi;
00102     sun_azimuth = sun_azimuth*invrad;
00103     // Parallax Correction
00104     dParallax=(dEarthMeanRadius/dAstronomicalUnit) *vcl_sin(sun_elevation);
00105     sun_elevation=(sun_elevation + dParallax)*invrad;
00106     sun_elevation = 90.0-sun_elevation;//angle from local horizon
00107   }
00108 }