ROPO
Loading...
Searching...
No Matches
fmi_sunpos.c
1
22#include <stdio.h>
23#include <stdlib.h>
24#include <string.h>
25#include <math.h>
26#include <time.h>
27#include <sys/types.h>
28#include <limits.h>
29#define KOR 0
30#define VAN 1
31#define ANJ 2
32#define IKA 3
33#define KUO 4
34#define UTA 5
35#define LUO 6
36
37
38/* 2001 (C) Harri.Hohti@fmi.fi */
39
40static double DEG_TO_RAD,RAD_TO_DEG;
41
42typedef struct { double RA;
43 double Dec;
44 } Equatorial;
45
46typedef struct { double A;
47 double h;
48 } Horizontal;
49
50typedef struct { double Lat;
51 double Lon;
52 } Geographical;
53
54/* Julian date */
55static double JDate(char *datestr);
56
57/* 0-meridiaanin astronomic time (tähtiaika) */
58static double Sidereal_GW(double JD);
59
60/* Equatorial coordinates of the Sun */
61static Equatorial Solar_coords(double JD);
62
63/* Horizontal coordinates of the Sun (no refraction) */
64static Horizontal Solar_pos(double Sd0,
65 Geographical Geo,
66 Equatorial Equ);
67
68
69void sunPosition(double latDeg, double lonDeg,char *timestamp,double *azimuth,double *elevation){
70
71 double JD, Sd0;
72 Equatorial Equ;
73 Horizontal Hor;
74 Geographical Geo;
75
76 Geo.Lat = latDeg * DEG_TO_RAD;
77 Geo.Lon = lonDeg * DEG_TO_RAD;
78
79 JD = JDate(timestamp);
80
81 /* T�htiaika radiaaneina auringon tuntikulman laskentaa varten */
82 Sd0 = Sidereal_GW(JD)*15.0*DEG_TO_RAD;
83
84 /* Lasketaan auringon ekvatoriaalikoordinaatit ajan funktiona */
85 Equ = Solar_coords(JD);
86
87 /* Lasketaan auringon horisontaalikoordinaatit paikan ja ajan funktiona */
88 Hor = Solar_pos(Sd0,Geo,Equ);
89
90 *azimuth = Hor.A;
91 *elevation = Hor.h;
92
93}
94
95/*/ 10.11.2010 */
96int mainOld(int argc, char *argv[])
97{
98
99 double JD,Sd0;
100 int r=0;
101 Equatorial Equ;
102 Horizontal Hor;
103 Geographical Geo;
104
105 /* Tutkien latitudit ja longitudit radiaaneina */
106 char radar[7][4]={"KOR","VAN","ANJ","IKA","KUO","UTA","LUO"};
107 Geographical Radars[7]={{1.0492337686,0.37757289484},
108 {1.0518517625,0.43400520732},
109 {1.0629055144,0.47298422728},
110 {1.0780317013,0.40258928078},
111 {1.0969394348,0.47792932682},
112 {1.1303915788,0.45931248147},
113 {1.1716977045,0.46949356877}};
114
115 DEG_TO_RAD=M_PI/180.0;
116 RAD_TO_DEG=180.0/M_PI;
117
118 /* JD=JDate(argv[1]); */
119
120 switch(argc){
121 case 4:
122 /* sunpos 200203011547 35.0 62.0 */
123 JD=JDate(argv[1]);
124 Geo.Lat=atof(argv[2])*DEG_TO_RAD;
125 Geo.Lon=atof(argv[3])*DEG_TO_RAD;
126 break;
127 case 2:
128 /* sunpos KOR200203011547 */
129 JD=JDate(&argv[1][3]);
130 for(r=0;r<7;r++) if(strncasecmp(argv[1],radar[r],3)==0) break;
131 Geo=Radars[r];
132 break;
133 case 3:
134 /* sunpos 200203011547 KOR */
135 JD=JDate(argv[1]);
136 for(r=0;r<7;r++) if(strncasecmp(argv[2],radar[r],3)==0) break;
137 Geo=Radars[r];
138 break;
139 default:
140 fprintf(stderr,"%s : illegal parameter count: %d\n",argv[0],argc);
141 fprintf(stderr,"usage:\n");
142 fprintf(stderr,"%s KOR200203011547\n",argv[0]);
143 fprintf(stderr,"%s 200203011547 KOR\n",argv[0]);
144 fprintf(stderr,"%s 200203011547 35.0 62.0 \n",argv[0]);
145 return (-1);
146 }
147
148 if(r==7) return(1);
149
150 /* T�htiaika radiaaneina auringon tuntikulman laskentaa varten */
151 Sd0=Sidereal_GW(JD)*15.0*DEG_TO_RAD;
152
153 /* Lasketaan auringon ekvatoriaalikoordinaatit ajan funktiona */
154 Equ=Solar_coords(JD);
155
156 /* Lasketaan auringon horisontaalikoordinaatit paikan ja ajan funktiona */
157 Hor=Solar_pos(Sd0,Geo,Equ);
158
159 /* printf("A = %.4f\nh = %.4f\n",Hor.A,Hor.h); */
160 printf("%.1f\t%.1f\n",Hor.A,Hor.h);
161
162 return(0);
163}
164
165static double Sidereal_GW(double JD)
166{
167 double JD0,d,T,SDR,SD_0GW,SD_GW;
168
169 JD0=(double)((long int)(JD-0.5));
170 d=(JD-0.5-JD0)*24.0;
171 T=(JD0+0.5-2415020.0)/36525.0;
172
173 /* printf("%f %.16f\n",JD0+0.5,T); */
174
175 /* Lasketaan t�htiaika 0-meridiaanilla 0 UTC */
176 SDR=0.276919398+100.0021359*T+0.000001075*T*T;
177 SD_0GW=(SDR-(double)((long int)SDR))*24.0;
178
179 /* ja lis�t��n siihen todellinen aika t�htiajaksi muunnettuna,
180 jotta saadaan t�htiaika halutulle ajanhetkelle */
181 SD_GW=SD_0GW+d*1.002737908;
182
183 return(SD_GW);
184}
185
186
187double JDate(char *datestr)
188{
189 int YY,MM,DD,hh,mm,A,B,y,m;
190 double d,JD=0.0,fy,fm;
191
192 sscanf(datestr,"%4d%2d%2d%2d%2d",&YY,&MM,&DD,&hh,&mm);
193
194 y=YY;
195 m=MM;
196 if(MM<3)
197 {
198 y=YY-1;
199 m=MM+12;
200 }
201 A=y/100;
202 B=2-A+A/4;
203 fy=(double)y;
204 fm=(double)m;
205
206 d=(double)DD+(double)hh/24.0+(double)mm/1440.0;
207 JD=(double)((long int)(fy*365.25))+(double)((long int)(30.6001*(fm+1)))
208 +d+1720994.5+(double)B;
209
210 return(JD);
211}
212
213static Equatorial Solar_coords(double JD)
214{
215
216 double T,TT,L,M,e,OL,C,W,y,x,OA,ep;
217 Equatorial Equ;
218
219 T=(JD-2415020.0)/36525.0;
220 TT=T*T;
221
222 /* Auringon keskilongitudi */
223 L=279.69668+36000.76892*T+0.0003025*TT;
224
225 /* Auringon keskianomalia */
226 M=358.47583+35999.04975*T-0.00015*TT-0.0000033*T*TT;
227
228 /* Maan radan eksentrisyys */
229 e=0.01675104-0.0000418*T-0.000000126*TT;
230
231 /* printf("M=%.12f\nL=%.12f\ne=%.12f\n",M,L,e); */
232
233 M*=DEG_TO_RAD;
234
235 /* Auringon "equation of the center" */
236 C=(1.919460-0.004789*T-0.000014*TT)*sin(M);
237 C+=(0.020094-0.0001*T)*sin(2*M);
238 C+=0.000293*sin(3*M);
239
240 /* Auringon longitudi */
241 OL=L+C;
242 /* printf("OL=%.12f\nC=%.12f\n",OL,C); */
243
244 /* Maan nousevan solmun longitudi */
245 W=(259.18-1934.142*T)*DEG_TO_RAD;
246
247 /* Auringon longitudi korjattuna todelliseen epookkiin */
248 OA=OL-0.00569-0.00479*sin(W);
249 /* printf("Oapp=%.12f\n",OA); */
250
251 OA*=DEG_TO_RAD;
252
253 /* ekliptikan kaltevuus */
254 ep=23.452294-0.0130125*T-0.00000164*TT+0.000000503*T*TT+0.00256*cos(W);
255 ep*=DEG_TO_RAD;
256
257 /* Muunnos auringon ekliptikaalisista koordinaateista
258 (longitudi lambda, latitudi beta (auringolle < 1.2"))
259 ekvatoriaalisiin, rektaskensioon ja deklinaatioon */
260
261 y=cos(ep)*sin(OA);
262 x=cos(OA);
263
264 Equ.Dec=asin(sin(ep)*sin(OA));
265 Equ.RA=atan2(y,x);
266
267 /* printf("RA=%.12f\nDec=%.12f\n",Equ.RA*RAD_TO_DEG/15.0,Equ.Dec*RAD_TO_DEG); */
268
269 return(Equ);
270}
271
272
273static Horizontal Solar_pos(double Sd0,Geographical Geo,
274 Equatorial Equ)
275{
276
277 double H,x,y;
278 Horizontal Hor;
279
280 /* Tuntikulma H on 0-meridiaanin t�htiaika + longitudi - rektaskensio */
281 H=Sd0+Geo.Lon-Equ.RA;
282
283 /* Muunnos ekvatoriaalikoordinaateista halutun latitudin ja
284 longitudin horisontaalisiin koordinaatteihin */
285
286 y=sin(H);
287 x=cos(H)*sin(Geo.Lat)-tan(Equ.Dec)*cos(Geo.Lat);
288
289 /* Lis�t��n 180, jotta atsimuutti alkaisi pohjoisesta */
290 Hor.A=atan2(y,x)*RAD_TO_DEG+180.0;
291
292 Hor.h=(asin(sin(Geo.Lat)*sin(Equ.Dec)
293 +cos(Geo.Lat)*cos(Equ.Dec)*cos(H)))*RAD_TO_DEG;
294
295 if(Hor.A<0.0) Hor.A+=360.0;
296
297 return(Hor);
298}
299
300
Generated by doxygen 1.11.0