/[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 2469 by jpye, Sat May 28 07:17:35 2011 UTC revision 2470 by jpye, Sat May 28 09:12:38 2011 UTC
# Line 1  Line 1 
1    /*  ASCEND modelling environment
2        Copyright (C) 2011 Carnegie Mellon University
3    
4        This program is free software; you can redistribute it and/or modify
5        it under the terms of the GNU General Public License as published by
6        the Free Software Foundation; either version 2, or (at your option)
7        any later version.
8    
9        This program is distributed in the hope that it will be useful,
10        but WITHOUT ANY WARRANTY; without even the implied warranty of
11        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12        GNU General Public License for more details.
13    
14        You should have received a copy of the GNU General Public License
15        along with this program; if not, write to the Free Software
16        Foundation, Inc., 59 Temple Place - Suite 330,
17        Boston, MA 02111-1307, USA.
18    *//** @FILE
19        This file's code is based on the source code given in the publication
20        R Grena (2008), An algorithm for the computation of the solar position,
21        Solar Energy (82), pp 462-470.
22    
23        The original code was in C++ and returned several intermediate results.
24        This modified version returns only the zenith and azimuth angles for
25        given date/time.
26    */
27    
28  #include "sunpos_grena.h"  #include "sunpos_grena.h"
29    #include <math.h>
30    
31  /* WARNING: this is a FIRST DRAFT of the code and HAS NOT BEEN CHECKED yet. */  /* we have converted this file from C++ to C. not so classy any more ;-) */
32    
33  void  
34  SunCoord::Calculate(){  double SunPos_calc_time(SunPos *S, double UT, int Day, int Month, int Year, double Delta_t){
35    
36      // calculation of JD and JDE      // calculation of JD and JDE
37      double dYear, dMonth;      double dYear, dMonth;
38      if(Month <= 2){      if(sunpos->Month <= 2){
39          dYear = (double)Year - 1.;          dYear = (double)S->Year - 1.;
40          dMonth = (double)Month + 12.;          dMonth = (double)S->Month + 12.;
41      }else{      }else{
42          dYear = (double)Year;          dYear = (double)S->Year;
43          dMonth = (double)Month;          dMonth = (double)S->Month;
44      }      }
45    
46      double JD_t = (double)trunc(365.25 * (dYear - 2000))      double JD_t = (double)trunc(365.25 * (dYear - 2000))
47          + (double)trunc(30.6001 * (dMonth + 1))          + (double)trunc(30.6001 * (dMonth + 1))
48          + (double)Day + UT/24. - 1158.5;          + (double)S->Day + S->UT/24. - 1158.5;
49    
50        double t = JD_t + sunpos->Delta_t/86400;
51    
52        S->t = t;
53    }
54    
55    
56    void SunPos_set_lat_long(SunPos *S, double latitude, double longitude){
57        S->latitude = latitude;
58        S->longitude = longitude;
59    }
60    
     double t = JD_t + Delta_t/86400;  
61    
62  #if 0  void SunPos_set_pressure_temp(SunPos *S, double p, double T){
63      // standard JD and JDE (useless for the computation, they are computed for      S->p = p;
64      // completeness:      S->T = T;
65    }
66    
67      // (omitted)  void SunPos_set_time(SunPos *S, double t){
68  #endif      S->t = t;
69    }
70    
71    
72    void SunPos_calc_zen_azi(SunPos *S, double *zenith, double *azimuth){
73        double t = S->t;
74        double HourAngle;
75        double TopocRightAscension;
76        double TopocDeclination;
77        double TopocHourAngle;
78        double Elevation_no_refrac;
79        double RefractionCorrection;
80    
81      // HELIOCENTRIC LONGITUDE      // HELIOCENTRIC LONGITUDE
82    
# Line 76  SunCoord::Calculate(){ Line 126  SunCoord::Calculate(){
126      // local hour angle of the sun      // local hour angle of the sun
127    
128      HourAngle = 6.30038809903 * JD_t + 4.8824623 + delta_psi * 0.9174      HourAngle = 6.30038809903 * JD_t + 4.8824623 + delta_psi * 0.9174
129          + ObserverLongitude - RightAscension;          + S->longitude - RightAscension;
130    
131      // to obtain the local hour angle in the range [0,2pi] uncomment:      // to obtain the local hour angle in the range [0,2pi] uncomment:
132      // HourAngle = fmod(HourAngle,2*PI);      // HourAngle = fmod(HourAngle,2*PI);
133    
134        double c_lat = cos(S->latitude);
135        double s_lat = sin(S->latitude);
136        double c_H = cos(HourAngle);
137        double s_H = sin(HourAngle);
138    
139      // parallax correction to right ascension      // parallax correction to right ascension
140    
141      double d_alpha = -4.26e-5 * c_lat * s_H;      double d_alpha = -4.26e-5 * c_lat * s_H;
# Line 111  SunCoord::Calculate(){ Line 166  SunCoord::Calculate(){
166    
167      // local coordinates of the sun      // local coordinates of the sun
168    
169      Zenith = PI/2 - Elevation_no_refrac - RefractionCorrection;      *zenith = PI/2 - Elevation_no_refrac - RefractionCorrection;
170    
171      Azimuth = atan2(s_H_corr, c_H_corr*s_lat - s_delta_corr/c_delta_corr*c_lat);      *azimuth = atan2(s_H_corr, c_H_corr*s_lat - s_delta_corr/c_delta_corr*c_lat);
172  }  }
173    
Legend:
Removed from v.2469  
changed lines
  Added in v.2470
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22