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

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