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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2469 - (show annotations) (download) (as text)
Sat May 28 07:17:35 2011 UTC (11 years, 1 month ago) by jpye
File MIME type: text/x-csrc
File size: 3288 byte(s)
First draft of Grena code, still needs to be checked.
1 #include "sunpos_grena.h"
2
3 /* WARNING: this is a FIRST DRAFT of the code and HAS NOT BEEN CHECKED yet. */
4
5 void
6 SunCoord::Calculate(){
7
8 // calculation of JD and JDE
9 double dYear, dMonth;
10 if(Month <= 2){
11 dYear = (double)Year - 1.;
12 dMonth = (double)Month + 12.;
13 }else{
14 dYear = (double)Year;
15 dMonth = (double)Month;
16 }
17
18 double JD_t = (double)trunc(365.25 * (dYear - 2000))
19 + (double)trunc(30.6001 * (dMonth + 1))
20 + (double)Day + UT/24. - 1158.5;
21
22 double t = JD_t + Delta_t/86400;
23
24 #if 0
25 // standard JD and JDE (useless for the computation, they are computed for
26 // completeness:
27
28 // (omitted)
29 #endif
30
31 // HELIOCENTRIC LONGITUDE
32
33 // linear increase + annual harmonic
34
35 double ang = 1.72019e-2 * t - 0.0563;
36 HeliocLongitude = 1.740940 + 1.7202768683e-2 * t + 3.34118e-2 * sin(ang) + 3.488e-4 * sin(2*ang);
37
38 // moon perturbation
39
40 HeliocLongitude += 3.13e-5 * sin(2.127730e-1*t - 0.585);
41
42 // harmonic correction
43
44 HeliocLongitude += 1.26e-5 * sin(4.243e-3 * t + 1.46)
45 + 2.35e-5 * sin(1.0727e-2 * t + 0.72)
46 + 2.76e-5 * sin(1.5799e-2 * t + 2.35)
47 + 2.75e-5 * sin(2.1551e-2 * t - 1.98)
48 + 1.26e-5 * sin(3.1490e-2 * t - 0.80);
49
50 // polynomial correction
51
52 double t2 = t/1000;
53 HeliocLongitude += ((( -2.30796e-7 * t2 + 3.7976e-6) * t2 - 2.0458e-5) * t + 3.976e-5) * t2*t2;
54
55 // to obtain obtain Heliocentric longitude in the range [0,2pi] uncomment:
56 // HeliocLongitude = fmod(HeliocLongitude, 2*PI);
57
58 // END HELIOCENTRIC LONGITUDE CALCULATION
59
60 // Correction to longitude due to nutation
61
62 double delta_psi = 8.33e-5 * sin(9.252e-4 * t - 1.173);
63
64 // Earth axis inclination
65
66 double epsilon = -6.21e-9 * t + 0.409086 + 4.46e-5 * sin(9.252e-4 * t + 0.397);
67
68 // Geocentric global solar coordinates:
69
70 GeocSolarLongitude = HeliocLongitude + PI + delta_psi - 9.932e-5;
71
72 double s_lambda = sin(GeocSolarLongitude);
73
74 RightAscension = atan2(s_lambda * cos(epsilon), cos(GeocSolarLongitude));
75
76 // local hour angle of the sun
77
78 HourAngle = 6.30038809903 * JD_t + 4.8824623 + delta_psi * 0.9174
79 + ObserverLongitude - RightAscension;
80
81 // to obtain the local hour angle in the range [0,2pi] uncomment:
82 // HourAngle = fmod(HourAngle,2*PI);
83
84 // parallax correction to right ascension
85
86 double d_alpha = -4.26e-5 * c_lat * s_H;
87 TopocRightAscension = RightAscension + d_alpha;
88 TopocHourAngle = HourAngle - d_alpha;
89
90 // parallax correction to declination:
91
92 TopocDeclination = Declination - 4.26e-5 * (s_lat - Declination * c_lat);
93 double s_delta_corr = sin(TopocDeclination);
94 double c_delta_corr = cos(TopocDeclination);
95 double c_H_corr = c_H + d_alpha * s_H;
96 double s_H_corr = s_H - d_alpha * c_H;
97
98 // solar elevation angle, without refraction correction
99 Elevation_no_refrac = asin(s_lat * s_delta_corr + c_lat * c_delta_corr * c_H_corr);
100
101 // refraction correction: it is calculated only
102 // if Elevation_no_refract > elev_min
103
104 const double elev_min = -0.01;
105
106 if(Elevation_no_refrac > elev_min){
107 RefractionCorrection = 0.084217 * Pressure / (273 + Temperature)
108 / tan(Elevation_no_refrac + 0.0031376 / (Elevation_no_refrac + 0.089186));
109 }else{
110 RefractionCorrection = 0;
111
112 // local coordinates of the sun
113
114 Zenith = PI/2 - Elevation_no_refrac - RefractionCorrection;
115
116 Azimuth = atan2(s_H_corr, c_H_corr*s_lat - s_delta_corr/c_delta_corr*c_lat);
117 }
118

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22