/[ascend]/trunk/models/johnpye/grena/sunpos_grena.c
ViewVC logotype

Diff of /trunk/models/johnpye/grena/sunpos_grena.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2470 by jpye, Sat May 28 09:12:38 2011 UTC revision 2471 by jpye, Sat May 28 10:24:10 2011 UTC
# Line 26  Line 26 
26  */  */
27    
28  #include "sunpos_grena.h"  #include "sunpos_grena.h"
 #include <math.h>  
29    
30  /* we have converted this file from C++ to C. not so classy any more ;-) */  /* we have converted this file from C++ to C. not so classy any more ;-) */
31    
32    void SunPos_calc_time(SunPos *S, double UT, int Day, int Month, int Year, double Delta_t){
 double SunPos_calc_time(SunPos *S, double UT, int Day, int Month, int Year, double Delta_t){  
33    
34      // calculation of JD and JDE      // calculation of JD and JDE
35      double dYear, dMonth;      double dYear, dMonth;
36      if(sunpos->Month <= 2){      if(Month <= 2){
37          dYear = (double)S->Year - 1.;          dYear = (double)Year - 1.;
38          dMonth = (double)S->Month + 12.;          dMonth = (double)Month + 12.;
39      }else{      }else{
40          dYear = (double)S->Year;          dYear = (double)Year;
41          dMonth = (double)S->Month;          dMonth = (double)Month;
42      }      }
43    
44      double JD_t = (double)trunc(365.25 * (dYear - 2000))      // universal time
45        S->t_G = (double)trunc(365.25 * (dYear - 2000))
46          + (double)trunc(30.6001 * (dMonth + 1))          + (double)trunc(30.6001 * (dMonth + 1))
47          + (double)S->Day + S->UT/24. - 1158.5;          + (double)Day + UT/24. - 1158.5;
48    
49        S->Delta_t = Delta_t;
50    }
51    
     double t = JD_t + sunpos->Delta_t/86400;  
52    
53      S->t = t;  void SunPos_set_time(SunPos *S, double t_G, double Delta_t){
54        S->t_G = t_G;
55        S->Delta_t = Delta_t;
56  }  }
57    
58    
# Line 64  void SunPos_set_pressure_temp(SunPos *S, Line 67  void SunPos_set_pressure_temp(SunPos *S,
67      S->T = T;      S->T = T;
68  }  }
69    
70  void SunPos_set_time(SunPos *S, double t){  /*
71      S->t = t;      Note: there is some confusion in nomenclature in the paper regarding
72  }      'JD_t' and 't_G' which are assumed to be the same time. This seems to John
73        to be correct.
74    */
75  void SunPos_calc_zen_azi(SunPos *S, double *zenith, double *azimuth){  void SunPos_calc_zen_azi(SunPos *S, double *zenith, double *azimuth){
     double t = S->t;  
76      double HourAngle;      double HourAngle;
77      double TopocRightAscension;      double TopocRightAscension;
78      double TopocDeclination;      double TopocDeclination;
# Line 80  void SunPos_calc_zen_azi(SunPos *S, doub Line 82  void SunPos_calc_zen_azi(SunPos *S, doub
82    
83      // HELIOCENTRIC LONGITUDE      // HELIOCENTRIC LONGITUDE
84    
85        double t = S->t_G + S->Delta_t / 86400;
86    
87      // linear increase + annual harmonic      // linear increase + annual harmonic
88    
89      double ang = 1.72019e-2 * t - 0.0563;      double ang = 1.72019e-2 * t - 0.0563;
90      HeliocLongitude = 1.740940 + 1.7202768683e-2 * t + 3.34118e-2 * sin(ang) + 3.488e-4 * sin(2*ang);      double HeliocLongitude = 1.740940 + 1.7202768683e-2 * t + 3.34118e-2 * sin(ang) + 3.488e-4 * sin(2*ang);
91    
92      // moon perturbation      // moon perturbation
93    
# Line 117  void SunPos_calc_zen_azi(SunPos *S, doub Line 121  void SunPos_calc_zen_azi(SunPos *S, doub
121    
122      // Geocentric global solar coordinates:      // Geocentric global solar coordinates:
123    
124      GeocSolarLongitude = HeliocLongitude + PI + delta_psi - 9.932e-5;      double GeocSolarLongitude = HeliocLongitude + PI + delta_psi - 9.932e-5;
125    
126      double s_lambda = sin(GeocSolarLongitude);      double s_lambda = sin(GeocSolarLongitude);
127    
128      RightAscension = atan2(s_lambda * cos(epsilon), cos(GeocSolarLongitude));      double RightAscension = atan2(s_lambda * cos(epsilon), cos(GeocSolarLongitude));
129    
130        double Declination = asin(sin(epsilon) * s_lambda);
131    
132      // local hour angle of the sun      // local hour angle of the sun
133    
134      HourAngle = 6.30038809903 * JD_t + 4.8824623 + delta_psi * 0.9174      HourAngle = 6.30038809903 * S->t_G + 4.8824623 + delta_psi * 0.9174
135          + S->longitude - RightAscension;          + S->longitude - RightAscension;
136    
137      // to obtain the local hour angle in the range [0,2pi] uncomment:      // to obtain the local hour angle in the range [0,2pi] uncomment:
# Line 159  void SunPos_calc_zen_azi(SunPos *S, doub Line 165  void SunPos_calc_zen_azi(SunPos *S, doub
165      const double elev_min = -0.01;      const double elev_min = -0.01;
166            
167      if(Elevation_no_refrac > elev_min){      if(Elevation_no_refrac > elev_min){
168          RefractionCorrection = 0.084217 * Pressure / (273 + Temperature)          RefractionCorrection = 0.084217 * S->p / (273 + S->T)
169              / tan(Elevation_no_refrac + 0.0031376 / (Elevation_no_refrac + 0.089186));              / tan(Elevation_no_refrac + 0.0031376 / (Elevation_no_refrac + 0.089186));
170      }else{      }else{
171          RefractionCorrection = 0;          RefractionCorrection = 0;
172        }
173    
174      // local coordinates of the sun      // local coordinates of the sun
175    
Legend:
Removed from v.2470  
changed lines
  Added in v.2471
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22