1 /* SWISSEPH
2
3 Heliacal risings and related calculations
4
5 Author: Victor Reijs
6 This program code is a translation of part of:
7 Victor Reijs' software ARCHAEOCOSMO (archaeoastronomy and
8 geodesy functions),
9 http://www.iol.ie/~geniet/eng/archaeocosmoprocedures.htm
10
11 Translation from VB into C by Dieter Koch
12
13 Problem reports can be sent to victor.reijs@gmail.com or dieter@astro.ch
14
15 Copyright (c) Victor Reijs, 2008
16 Copyright (C) 1997 - 2021 Astrodienst AG, Switzerland. All rights reserved.
17
18 License conditions
19 ------------------
20
21 This file is part of Swiss Ephemeris.
22
23 Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND. No author
24 or distributor accepts any responsibility for the consequences of using it,
25 or for whether it serves any particular purpose or works at all, unless he
26 or she says so in writing.
27
28 Swiss Ephemeris is made available by its authors under a dual licensing
29 system. The software developer, who uses any part of Swiss Ephemeris
30 in his or her software, must choose between one of the two license models,
31 which are
32 a) GNU Affero General Public License (AGPL)
33 b) Swiss Ephemeris Professional License
34
35 The choice must be made before the software developer distributes software
36 containing parts of Swiss Ephemeris to others, and before any public
37 service using the developed software is activated.
38
39 If the developer choses the AGPL software license, he or she must fulfill
40 the conditions of that license, which includes the obligation to place his
41 or her whole software project under the AGPL or a compatible license.
42 See https://www.gnu.org/licenses/agpl-3.0.html
43
44 If the developer choses the Swiss Ephemeris Professional license,
45 he must follow the instructions as found in http://www.astro.com/swisseph/
46 and purchase the Swiss Ephemeris Professional Edition from Astrodienst
47 and sign the corresponding license contract.
48
49 The License grants you the right to use, copy, modify and redistribute
50 Swiss Ephemeris, but only under certain conditions described in the License.
51 Among other things, the License requires that the copyright notices and
52 this notice be preserved on all copies.
53
54 Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl
55
56 The authors of Swiss Ephemeris have no control or influence over any of
57 the derived works, i.e. over software or services created by other
58 programmers which use Swiss Ephemeris functions.
59
60 The names of the authors or of the copyright holder (Astrodienst) must not
61 be used for promoting any software, product or service which uses or contains
62 the Swiss Ephemeris. This copyright notice is the ONLY place where the
63 names of the authors can legally appear, except in cases where they have
64 given special permission in writing.
65
66 The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
67 for promoting such software, products or services.
68 */
69
70
71 #include "swephexp.h"
72 #include "sweph.h"
73 #include "swephlib.h"
74 #include <sys/stat.h>
75
76 #define PLSV 0 /*if Planet, Lunar and Stellar Visibility formula is needed PLSV=1*/
77 #define criticalangle 0.0 /*[deg]*/
78 #define BNIGHT 1479.0 /*[nL]*/
79 #define BNIGHT_FACTOR 1.0
80 #define PI M_PI
81 #define Min2Deg (1.0 / 60.0)
82 #define SWEHEL_DEBUG 0
83 #define DONE 1
84 #define MaxTryHours 4
85 #define TimeStepDefault 1
86 #define LocalMinStep 8
87
88 /* time constants */
89 #define Y2D 365.25 /*[Day]*/
90 #define D2Y (1 / Y2D) /*[Year]*/
91 #define D2H 24.0 /*[Hour]*/
92 #define H2S 3600.0 /*[sec]*/
93 #define D2S (D2H * H2S) /*[sec]*/
94 #define S2H (1.0 / H2S) /*[Hour]*/
95 #define JC2D 36525.0 /*[Day]*/
96 #define M2S 60.0 /*[sec]*/
97
98 /* Determines which algorimths are used*/
99 #define REFR_SINCLAIR 0
100 #define REFR_BENNETTH 1
101 #define FormAstroRefrac REFR_SINCLAIR /*for Astronomical refraction can be "bennetth" or "sinclair"*/
102 #define GravitySource 2 /*0=RGO, 1=Wikipedia,2=Exp. Suppl. 1992,3=van der Werf*/
103 #define REarthSource 1 /*0=RGO (constant), 1=WGS84 method*/
104
105 #define StartYear 1820 /*[year]*/
106 #define Average 1.80546834626888 /*[msec/cy]*/
107 #define Periodicy 1443.67123144531 /*[year]*/
108 #define Amplitude 3.75606495492684 /*[msec]*/
109 #define phase 0 /*[deg]*/
110 #define MAX_COUNT_SYNPER 5 /* search within 10 synodic periods */
111 #define MAX_COUNT_SYNPER_MAX 1000000 /* high, so there is not max count */
112 #define AvgRadiusMoon (15.541 / 60) /* '[Deg] at 2007 CE or BCE*/
113
114 /* WGS84 ellipsoid constants
115 * http://w3sli.wcape.gov.za/Surveys/Mapping/wgs84.htm*/
116 #define Ra 6378136.6 /*'[m]*/
117 #define Rb 6356752.314 /*'[m]*/
118
119 /* choices in Schaefer's model */
120 #define nL2erg (1.02E-15)
121 #define erg2nL (1 / nL2erg) /*erg2nL to nLambert*/
122 #define MoonDistance 384410.4978 /*[km]*/
123 #define scaleHwater 3000.0 /*[m] Ricchiazzi [1997] 8200 Schaefer [2000]*/
124 #define scaleHrayleigh 8515.0 /*[m] Su [2003] 8200 Schaefer [2000]*/
125 #define scaleHaerosol 3745.0 /*m Su [2003] 1500 Schaefer [2000]*/
126 #define scaleHozone 20000.0 /*[m] Schaefer [2000]*/
127 #define astr2tau 0.921034037197618 /*LN(10 ^ 0.4)*/
128 #define tau2astr 1 / astr2tau
129
130 /* meteorological constants*/
131 #define C2K 273.15 /*[K]*/
132 #define DELTA 18.36
133 #define TempNulDiff 0.000001
134 #define PressRef 1000 /*[mbar]*/
135 #define MD 28.964 /*[kg] Mol weight of dry air van der Werf*/
136 #define MW 18.016 /*[kg] Mol weight of water vapor*/
137 #define GCR 8314.472 /*[L/kmol/K] van der Werf*/
138 #define LapseSA 0.0065 /*[K/m] standard atmosphere*/
139 #define LapseDA 0.0098 /*[K/m] dry adiabatic*/
140
141 /* lowest apparent altitude to provide*/
142 #define LowestAppAlt -3.5 /*[Deg]*/
143
144 /*optimization delta*/
145 #define epsilon 0.001
146 /* for Airmass usage*/
147 #define staticAirmass 0 /* use staticAirmass=1 instead depending on difference k's*/
148
149 /* optic stuff */
150 #define GOpticMag 1 /*telescope magnification*/
151 #define GOpticTrans 0.8 /*telescope transmission*/
152 #define GBinocular 1 /*1-binocular 0=monocular*/
153 #define GOpticDia 50 /*telescope diameter [mm]*/
154
mymin(double a,double b)155 static double mymin(double a, double b)
156 {
157 if (a <= b)
158 return a;
159 return b;
160 }
161
mymax(double a,double b)162 static double mymax(double a, double b)
163 {
164 if (a >= b)
165 return a;
166 return b;
167 }
168
169 /*###################################################################*/
Tanh(double x)170 static double Tanh(double x)
171 {
172 return (exp(x) - exp(-x)) / (exp(x) + exp(-x));
173 }
174
175 /*
176 ' B [nL]
177 ' SN [-]
178 ' CVA [deg]
179 */
CVA(double B,double SN,int32 helflag)180 static double CVA(double B, double SN, int32 helflag)
181 {
182 /*Schaefer, Astronomy and the limits of vision, Archaeoastronomy, 1993*/
183 AS_BOOL is_scotopic = FALSE;
184 //if (B < BNIGHT)
185 if (B < 1394) /* use this value for BNIGHT to make the function continous */
186 is_scotopic = TRUE;
187 if (helflag & SE_HELFLAG_VISLIM_PHOTOPIC)
188 is_scotopic = FALSE;
189 if (helflag & SE_HELFLAG_VISLIM_SCOTOPIC)
190 is_scotopic = TRUE;
191 if (is_scotopic)
192 return mymin(900, 380 / SN * pow(10, (0.3 * pow(B, (-0.29))))) / 60.0 / 60.0;
193 else
194 return (40.0 / SN) * pow(10, (8.28 * pow(B, (-0.29)))) / 60.0 / 60.0;
195 }
196
197 /*
198 ' age [year]
199 ' B [nL]
200 ' PupilDia [mm]
201 */
PupilDia(double Age,double B)202 static double PupilDia(double Age, double B)
203 {
204 /* age dependancy from Garstang [2000]*/
205 return (0.534 - 0.00211 * Age - (0.236 - 0.00127 * Age) * Tanh(0.4 * log(B) / log(10) - 2.2)) * 10;
206 }
207
208 /*
209 'Input
210 ' Bback [nL]
211 ' kX [-]
212 ' Binocular [-]
213 ' OpticMag [-]
214 ' OpticDia [mm]
215 ' OpticTrans [-]
216 ' JDNDaysUT [JDN]
217 ' Age [Year]
218 ' SN [-]
219 ' ObjectName
220 ' TypeFactor [0=itensity factor 1=background factor]
221 'Output
222 ' OpticFactor [-]
223 */
OpticFactor(double Bback,double kX,double * dobs,double JDNDaysUT,char * ObjectName,int TypeFactor,int helflag)224 static double OpticFactor(double Bback, double kX, double *dobs, double JDNDaysUT, char *ObjectName, int TypeFactor, int helflag)
225 {
226 double Pst, CIb, CIi, ObjectSize, Fb, Fe, Fsc, Fci, Fcb, Ft, Fp, Fa, Fr, Fm;
227 double Age = dobs[0];
228 double SN = dobs[1], SNi;
229 double Binocular = dobs[2];
230 double OpticMag = dobs[3];
231 double OpticDia = dobs[4];
232 double OpticTrans = dobs[5];
233 AS_BOOL is_scotopic = FALSE;
234 JDNDaysUT += 0.0; /* currently not used, statement prevents compiler warning */
235 SNi = SN;
236 if (SNi <= 0.00000001) SNi = 0.00000001;
237 /* 23 jaar as standard from Garstang*/
238 Pst = PupilDia(23, Bback);
239 if (OpticMag == 1) { /*OpticMagn=1 means using eye*/
240 OpticTrans = 1;
241 OpticDia = Pst;
242 }
243 #if 0 /*is done in default_heliacal_parameters()*/
244 if (OpticMag == 0) { /*OpticMagn=0 (undefined) using eye*/
245 OpticTrans = 1;
246 OpticDia = Pst;
247 Binocular = 1;
248 OpticMag = 1;
249 }
250 #endif
251 /* Schaefer, Astronomy and the limits of vision, Archaeoastronomy, 1993*/
252 CIb = 0.7; /* color of background (from Ben Sugerman)*/
253 CIi = 0.5; /* Color index for white (from Ben Sugerman), should be function of ObjectName*/
254 ObjectSize = 0;
255 if (strcmp(ObjectName, "moon") == 0) {
256 /*ObjectSize and CI needs to be determined (depending on JDNDaysUT)*/
257 ;
258 }
259 Fb = 1;
260 if (Binocular == 0) Fb = 1.41;
261 //if (Bback < BNIGHT)
262 if (Bback < 1645) /* use this value for BNIGHT to make the function continuous */
263 is_scotopic = TRUE;
264 if (helflag & SE_HELFLAG_VISLIM_PHOTOPIC)
265 is_scotopic = FALSE;
266 if (helflag & SE_HELFLAG_VISLIM_SCOTOPIC)
267 is_scotopic = TRUE;
268 if (is_scotopic) {
269 Fe = pow(10, (0.48 * kX));
270 Fsc = mymin(1, (1 - pow(Pst / 124.4, 4)) / (1 - pow((OpticDia / OpticMag / 124.4), 4)));
271 Fci = pow(10, (-0.4 * (1 - CIi / 2.0)));
272 Fcb = pow(10, (-0.4 * (1 - CIb / 2.0)));
273 } else {
274 Fe = pow(10, (0.4 * kX));
275 Fsc = mymin(1, pow((OpticDia / OpticMag / Pst), 2) * (1 - exp(-pow((Pst / 6.2), 2))) / (1 - exp(-pow((OpticDia / OpticMag / 6.2), 2))));
276 Fci = 1;
277 Fcb = 1;
278 }
279 Ft = 1 / OpticTrans;
280 Fp = mymax(1, pow((Pst / (OpticMag * PupilDia(Age, Bback))), 2));
281 Fa = pow((Pst / OpticDia), 2);
282 Fr = (1 + 0.03 * pow((OpticMag * ObjectSize / CVA(Bback, SNi, helflag)), 2)) / pow(SNi, 2);
283 Fm = pow(OpticMag, 2);
284 #if SWEHEL_DEBUG
285 fprintf(stderr, "Pst=%f\n", Pst);
286 fprintf(stderr, "Fb =%f\n", Fb);
287 fprintf(stderr, "Fe =%f\n", Fe);
288 fprintf(stderr, "Ft =%f\n", Ft);
289 fprintf(stderr, "Fp =%f\n", Fp);
290 fprintf(stderr, "Fa =%f\n", Fa);
291 fprintf(stderr, "Fm =%f\n", Fm);
292 fprintf(stderr, "Fsc=%f\n", Fsc);
293 fprintf(stderr, "Fci=%f\n", Fci);
294 fprintf(stderr, "Fcb=%f\n", Fcb);
295 fprintf(stderr, "Fr =%f\n", Fr );
296 #endif
297 if (TypeFactor == 0)
298 return Fb * Fe * Ft * Fp * Fa * Fr * Fsc * Fci;
299 else
300 return Fb * Ft * Fp * Fa * Fm * Fsc * Fcb;
301 }
302
303 /*###################################################################
304 */
DeterObject(char * ObjectName)305 static int32 DeterObject(char *ObjectName)
306 {
307 char s[AS_MAXCH];
308 char *sp;
309 int32 ipl;
310 strcpy(s, ObjectName);
311 for (sp = s; *sp != '\0'; sp++)
312 *sp = tolower(*sp);
313 if (strncmp(s, "sun", 3) == 0)
314 return SE_SUN;
315 if (strncmp(s, "venus", 5) == 0)
316 return SE_VENUS;
317 if (strncmp(s, "mars", 4) == 0)
318 return SE_MARS;
319 if (strncmp(s, "mercur", 6) == 0)
320 return SE_MERCURY;
321 if (strncmp(s, "jupiter", 7) == 0)
322 return SE_JUPITER;
323 if (strncmp(s, "saturn", 6) == 0)
324 return SE_SATURN;
325 if (strncmp(s, "uranus", 6) == 0)
326 return SE_URANUS;
327 if (strncmp(s, "neptun", 6) == 0)
328 return SE_NEPTUNE;
329 if (strncmp(s, "moon", 4) == 0)
330 return SE_MOON;
331 if ((ipl = atoi(s)) > 0) {
332 ipl += SE_AST_OFFSET;
333 return ipl;
334 }
335 return -1;
336 }
337
338 #if 0
339 int32 call_swe_calc(double tjd, int32 ipl, int32 iflag, double *x, char *serr)
340 {
341 int32 retval = OK, ipli, i;
342 double dtjd;
343 static TLS double tjdsv[3];
344 static TLS double xsv[3][6];
345 static TLS int32 iflagsv[3];
346 ipli = ipl;
347 if (ipli > SE_MOON)
348 ipli = 2;
349 dtjd = tjd - tjdsv[ipli];
350 if (tjdsv[ipli] != 0 && iflag == iflagsv[ipli] && fabs(dtjd) < 5.0 / 1440.0) {
351 for (i = 0; i < 3; i++)
352 x[i] = xsv[ipli][i] + dtjd * xsv[ipli][i+3];
353 for (i = 3; i < 6; i++)
354 x[i] = xsv[ipli][i];
355 } else {
356 retval = swe_calc(tjd, ipl, iflag, x, serr);
357 tjdsv[ipli] = tjd;
358 iflagsv[ipli] = iflag;
359 for (i = 0; i < 6; i++)
360 xsv[ipli][i] = x[i];
361 }
362 return retval;
363 }
364 #endif
365
366 /* avoids problems with star name string that may be overwritten by
367 swe_fixstar() */
call_swe_fixstar(char * star,double tjd,int32 iflag,double * xx,char * serr)368 static int32 call_swe_fixstar(char *star, double tjd, int32 iflag, double *xx, char *serr)
369 {
370 int32 retval;
371 char star2[AS_MAXCH];
372 strcpy(star2, star);
373 retval = swe_fixstar(star2, tjd, iflag, xx, serr);
374 return retval;
375 }
376
377 /* avoids problems with star name string that may be overwritten by
378 swe_fixstar_mag() */
call_swe_fixstar_mag(char * star,double * mag,char * serr)379 static int32 call_swe_fixstar_mag(char *star, double *mag, char *serr)
380 {
381 int32 retval;
382 char star2[AS_MAXCH];
383 static TLS double dmag;
384 static TLS char star_save[AS_MAXCH];
385 if (strcmp(star, star_save) == 0) {
386 *mag = dmag;
387 return OK;
388 }
389 strcpy(star_save, star);
390 strcpy(star2, star);
391 retval = swe_fixstar_mag(star2, &dmag, serr);
392 *mag = dmag;
393 return retval;
394 }
395
396 /* avoids problems with star name string that may be overwritten by
397 swe_fixstar() */
call_swe_rise_trans(double tjd,int32 ipl,char * star,int32 helflag,int32 eventtype,double * dgeo,double atpress,double attemp,double * tret,char * serr)398 static int32 call_swe_rise_trans(double tjd, int32 ipl, char *star, int32 helflag, int32 eventtype, double *dgeo, double atpress, double attemp, double *tret, char *serr)
399 {
400 int32 retval;
401 int32 iflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
402 char star2[AS_MAXCH];
403 strcpy(star2, star);
404 retval = swe_rise_trans(tjd, ipl, star2, iflag, eventtype, dgeo, atpress, attemp, tret, serr);
405 return retval;
406 }
407
408 /*
409 * Written by Dieter Koch:
410 * Fast function for risings and settings of planets, can be used instead of
411 * swe_rise_trans(), which is much slower.
412 * For circumpolar and near-circumpolar planets use swe_rise_trans(), or
413 * generally use it for geographical latitudes higher than 58N/S.
414 * For fixed stars, swe_rise_trans() is fast enough.
415 */
calc_rise_and_set(double tjd_start,int32 ipl,double * dgeo,double * datm,int32 eventflag,int32 helflag,double * trise,char * serr)416 static int32 calc_rise_and_set(double tjd_start, int32 ipl, double *dgeo, double *datm, int32 eventflag, int32 helflag, double *trise, char *serr)
417 {
418 int retc = OK, i;
419 double sda, xs[6], xx[6], xaz[6], xaz2[6], dfac = 1/365.25;
420 double rdi, rh;
421 double tjd0 = tjd_start, tjdrise;
422 double tjdnoon = (int) tjd0 - dgeo[0] / 15.0 / 24.0;
423 int32 iflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
424 int32 epheflag = iflag;
425 iflag |= SEFLG_EQUATORIAL;
426 if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
427 iflag |= SEFLG_NONUT|SEFLG_TRUEPOS;
428 if (swe_calc_ut(tjd0, SE_SUN, iflag, xs, serr) == 0) {
429 if (serr != NULL)
430 strcpy(serr, "error in calc_rise_and_set(): calc(sun) failed ");
431 return ERR;
432 }
433 if (swe_calc_ut(tjd0, ipl, iflag, xx, serr) == 0) {
434 if (serr != NULL)
435 strcpy(serr, "error in calc_rise_and_set(): calc(sun) failed ");
436 return ERR;
437 }
438 tjdnoon -= swe_degnorm(xs[0] - xx[0])/360.0 + 0;
439 /* is planet above horizon or below? */
440 swe_azalt(tjd0, SE_EQU2HOR, dgeo, datm[0], datm[1], xx, xaz);
441 if (eventflag & SE_CALC_RISE) {
442 if (xaz[2] > 0) {
443 while (tjdnoon - tjd0 < 0.5) {/*printf("e");*/tjdnoon += 1;}
444 while (tjdnoon - tjd0 > 1.5) {/*printf("f");*/tjdnoon -= 1;}
445 } else {
446 while (tjdnoon - tjd0 < 0.0) {/*printf("g");*/tjdnoon += 1;}
447 while (tjdnoon - tjd0 > 1.0) {/*printf("h");*/tjdnoon -= 1;}
448 }
449 } else {
450 if (xaz[2] > 0) {
451 while (tjd0 - tjdnoon > 0.5) {/*printf("a");*/ tjdnoon += 1;}
452 while (tjd0 - tjdnoon < -0.5) {/*printf("b");*/ tjdnoon -= 1;}
453 } else {
454 while (tjd0 - tjdnoon > 0.0) {/*printf("c");*/ tjdnoon += 1;}
455 while (tjd0 - tjdnoon < -1.0) {/*printf("d");*/ tjdnoon -= 1;}
456 }
457 }
458 /* position of planet */
459 if (swe_calc_ut(tjdnoon, ipl, iflag, xx, serr) == ERR) {
460 if (serr != NULL)
461 strcpy(serr, "error in calc_rise_and_set(): calc(sun) failed ");
462 return ERR;
463 }
464 /* apparent radius of solar disk (ignoring refraction) */
465 rdi = 0;
466 if (ipl == SE_SUN)
467 rdi = asin(696000000.0 / 1.49597870691e+11 / xx[2]) / DEGTORAD;
468 else if (ipl == SE_MOON)
469 rdi = asin(1737000.0 / 1.49597870691e+11 / xx[2]) / DEGTORAD;
470 if (eventflag & SE_BIT_DISC_CENTER)
471 rdi = 0;
472 /* true altitude of sun, when it appears at the horizon */
473 /* refraction for a body visible at the horizon at 0m above sea,
474 * atmospheric temperature 10 deg C, atmospheric pressure 1013.25 is 34.5 arcmin*/
475 rh = -(34.5 / 60.0 + rdi);
476 /* semidiurnal arc of sun */
477 sda = acos(-tan(dgeo[1] * DEGTORAD) * tan(xx[1] * DEGTORAD)) * RADTODEG;
478 /* rough rising and setting times */
479 if (eventflag & SE_CALC_RISE)
480 tjdrise = tjdnoon - sda / 360.0;
481 else
482 tjdrise = tjdnoon + sda / 360.0;
483 /*ph->tset = tjd_start + sda / 360.0;*/
484 /* now calculate more accurate rising and setting times.
485 * use vertical speed in order to determine crossing of the horizon
486 * refraction of 34' and solar disk diameter of 16' = 50' = 0.84 deg */
487 iflag = epheflag|SEFLG_SPEED|SEFLG_EQUATORIAL;
488 if (ipl == SE_MOON)
489 iflag |= SEFLG_TOPOCTR;
490 if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
491 iflag |= SEFLG_NONUT|SEFLG_TRUEPOS;
492 for (i = 0; i < 2; i++) {
493 if (swe_calc_ut(tjdrise, ipl, iflag, xx, serr) == ERR) {
494 /*fprintf(stderr, "hev4 tjd=%f, ipl=%d, iflag=%d\n", tjdrise, ipl, iflag);*/
495 return ERR;
496 }
497 swe_azalt(tjdrise, SE_EQU2HOR, dgeo, datm[0], datm[1], xx, xaz);
498 xx[0] -= xx[3] * dfac;
499 xx[1] -= xx[4] * dfac;
500 swe_azalt(tjdrise - dfac, SE_EQU2HOR, dgeo, datm[0], datm[1], xx, xaz2);
501 tjdrise -= (xaz[1] - rh) / (xaz[1] - xaz2[1]) * dfac;
502 /*fprintf(stderr, "%f\n", ph->trise);*/
503 }
504 *trise = tjdrise;
505 return retc;
506 }
507
my_rise_trans(double tjd,int32 ipl,char * starname,int32 eventtype,int32 helflag,double * dgeo,double * datm,double * tret,char * serr)508 static int32 my_rise_trans(double tjd, int32 ipl, char* starname, int32 eventtype, int32 helflag, double *dgeo, double *datm, double *tret, char *serr)
509 {
510 int retc = OK;
511 if (starname != NULL && *starname != '\0')
512 ipl = DeterObject(starname);
513 /* for non-circumpolar planets we can use a faster algorithm */
514 /*if (!(helflag & SE_HELFLAG_HIGH_PRECISION) && ipl != -1 && fabs(dgeo[1]) < 58) {*/
515 if (ipl != -1 && fabs(dgeo[1]) < 63) {
516 retc = calc_rise_and_set(tjd, ipl, dgeo, datm, eventtype, helflag, tret, serr);
517 /* for stars and circumpolar planets we use a rigorous algorithm */
518 } else {
519 retc = call_swe_rise_trans(tjd, ipl, starname, helflag, eventtype, dgeo, datm[0], datm[1], tret, serr);
520 }
521 /* printf("%f, %f\n", tjd, *tret);*/
522 return retc;
523 }
524
525 /*###################################################################
526 ' JDNDaysUT [Days]
527 ' dgeo [array: longitude, latitude, eye height above sea m]
528 ' TempE [C]
529 ' PresE [mbar]
530 ' ObjectName (string)
531 ' RSEvent (1=rise, 2=set,3=up transit,4=down transit)
532 ' Rim [0=center,1=top]
533 ' RiseSet [Day]
534 */
RiseSet(double JDNDaysUT,double * dgeo,double * datm,char * ObjectName,int32 RSEvent,int32 helflag,int32 Rim,double * tret,char * serr)535 static int32 RiseSet(double JDNDaysUT, double *dgeo, double *datm, char *ObjectName, int32 RSEvent, int32 helflag, int32 Rim, double *tret, char *serr)
536 {
537 int32 eventtype = RSEvent, Planet, retval;
538 if (Rim == 0)
539 eventtype |= SE_BIT_DISC_CENTER;
540 Planet = DeterObject(ObjectName);
541 if (Planet != -1)
542 retval = my_rise_trans(JDNDaysUT, Planet, "", eventtype, helflag, dgeo, datm, tret, serr);
543 else
544 retval = my_rise_trans(JDNDaysUT, -1, ObjectName, eventtype, helflag, dgeo, datm, tret, serr);
545 return retval;
546 }
547
548 /*###################################################################
549 ' JDNDaysUT [Days]
550 ' actual [0= approximation, 1=actual]
551 ' SunRA [deg]
552 */
SunRA(double JDNDaysUT,int32 helflag,char * serr)553 static double SunRA(double JDNDaysUT, int32 helflag, char *serr)
554 {
555 int imon, iday, iyar, calflag = SE_GREG_CAL;
556 double dut;
557 static TLS double tjdlast;
558 static TLS double ralast;
559 helflag += 0; /* statement prevents compiler warning */
560 *serr = '\0';
561 if (JDNDaysUT == tjdlast)
562 return ralast;
563 #ifndef SIMULATE_VICTORVB
564 if (1) { /*helflag & SE_HELFLAG_HIGH_PRECISION) {*/
565 double tjd_tt;
566 double x[6];
567 int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
568 int32 iflag = epheflag | SEFLG_EQUATORIAL;
569 iflag |= SEFLG_NONUT | SEFLG_TRUEPOS;
570 tjd_tt = JDNDaysUT + swe_deltat_ex(JDNDaysUT, epheflag, serr);
571 if (swe_calc(tjd_tt, SE_SUN, iflag, x, serr) != ERR) {
572 ralast = x[0];
573 tjdlast = JDNDaysUT;
574 return ralast;
575 }
576 }
577 #endif
578 swe_revjul(JDNDaysUT, calflag, &iyar, &imon, &iday, &dut); /* this seems to be much faster than calling swe_revjul() ! Note: only because SunRA is called 1000s of times */
579 tjdlast = JDNDaysUT;
580 ralast = swe_degnorm((imon + (iday - 1) / 30.4 - 3.69) * 30);
581 /*ralast = (DatefromJDut(JDNDaysUT, 2) + (DatefromJDut(JDNDaysUT, 3) - 1) / 30.4 - 3.69) * 30;*/
582 return ralast;
583 }
584
585 /*###################################################################
586 ' Temp [C]
587 ' Kelvin [K]
588 */
Kelvin(double Temp)589 static double Kelvin(double Temp)
590 {
591 /*' http://en.wikipedia.org/wiki/Kelvin*/
592 return Temp + C2K;
593 }
594
595 /*###################################################################
596 ' AppAlt [deg]
597 ' TempE [C]
598 ' PresE [mbar]
599 ' TopoAltitudefromAppAlt [deg]
600 */
TopoAltfromAppAlt(double AppAlt,double TempE,double PresE)601 static double TopoAltfromAppAlt(double AppAlt, double TempE, double PresE)
602 {
603 double R = 0;
604 double retalt = 0;
605 if (AppAlt >= LowestAppAlt) {
606 if (AppAlt > 17.904104638432)
607 R = 0.97 / tan(AppAlt * DEGTORAD);
608 else
609 R = (34.46 + 4.23 * AppAlt + 0.004 * AppAlt * AppAlt) / (1 + 0.505 * AppAlt + 0.0845 * AppAlt * AppAlt);
610 R = (PresE - 80) / 930 / (1 + 0.00008 * (R + 39) * (TempE - 10)) * R;
611 retalt = AppAlt - R * Min2Deg;
612 } else {
613 retalt = AppAlt;
614 }
615 return retalt;
616 }
617
618 /*###################################################################
619 ' TopoAlt [deg]
620 ' TempE [C]
621 ' PresE [mbar]
622 ' AppAltfromTopoAlt [deg]
623 ' call this instead of swe_azalt(), because it is faster (lower precision
624 ' is required)
625 */
AppAltfromTopoAlt(double TopoAlt,double TempE,double PresE,int32 helflag)626 static double AppAltfromTopoAlt(double TopoAlt, double TempE, double PresE, int32 helflag)
627 {
628 /* using methodology of Newtown derivatives (analogue to what Swiss Emphemeris uses)*/
629 int i, nloop = 2;
630 double newAppAlt = TopoAlt;
631 double newTopoAlt = 0.0;
632 double oudAppAlt = newAppAlt;
633 double oudTopoAlt = newTopoAlt;
634 double verschil, retalt;
635 if (helflag & SE_HELFLAG_HIGH_PRECISION)
636 nloop = 5;
637 for (i = 0; i <= nloop; i++) {
638 newTopoAlt = newAppAlt - TopoAltfromAppAlt(newAppAlt, TempE, PresE);
639 /*newTopoAlt = newAppAlt - swe_refrac(newAppAlt, PresE, TempE, SE_CALC_APP_TO_TRUE);*/
640 verschil = newAppAlt - oudAppAlt;
641 oudAppAlt = newTopoAlt - oudTopoAlt - verschil;
642 if ((verschil != 0) && (oudAppAlt != 0))
643 verschil = newAppAlt - verschil * (TopoAlt + newTopoAlt - newAppAlt) / oudAppAlt;
644 else
645 verschil = TopoAlt + newTopoAlt;
646 oudAppAlt = newAppAlt;
647 oudTopoAlt = newTopoAlt;
648 newAppAlt = verschil;
649 }
650 retalt = TopoAlt + newTopoAlt;
651 if (retalt < LowestAppAlt)
652 retalt = TopoAlt;
653 return retalt;
654 }
655
656 /*###################################################################
657 ' TopoAlt [deg]
658 ' TopoDecl [deg]
659 ' Lat [deg]
660 ' HourAngle [hour]
661 */
HourAngle(double TopoAlt,double TopoDecl,double Lat)662 static double HourAngle(double TopoAlt, double TopoDecl, double Lat)
663 {
664 double Alti = TopoAlt * DEGTORAD;
665 double decli = TopoDecl * DEGTORAD;
666 double Lati = Lat * DEGTORAD;
667 double ha = (sin(Alti) - sin(Lati) * sin(decli)) / cos(Lati) / cos(decli);
668 if (ha < -1) ha = -1;
669 if (ha > 1) ha = 1;
670 /* from http://star-www.st-and.ac.uk/~fv/webnotes/chapt12.htm*/
671 return acos(ha) / DEGTORAD / 15.0;
672 }
673
674 /*###################################################################
675 ' JDNDaysUT [Days]
676 ' dgeo [array: longitude, latitude, eye height above sea m]
677 ' TempE [C]
678 ' PresE [mbar]
679 ' ObjectName [-]
680 ' Angle (0 = TopoAlt, 1 = Azi, 2=Topo Declination, 3=Topo Rectascension, 4=AppAlt,5=Geo Declination, 6=Geo Rectascension)
681 ' ObjectLoc [deg]
682 */
ObjectLoc(double JDNDaysUT,double * dgeo,double * datm,char * ObjectName,int32 Angle,int32 helflag,double * dret,char * serr)683 static int32 ObjectLoc(double JDNDaysUT, double *dgeo, double *datm, char *ObjectName, int32 Angle, int32 helflag, double *dret, char *serr)
684 {
685 double x[6], xin[3], xaz[3], tjd_tt;
686 int32 Planet;
687 int32 epheflag;
688 int32 iflag = SEFLG_EQUATORIAL;
689 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
690 iflag |= epheflag;
691 if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
692 iflag |= SEFLG_NONUT | SEFLG_TRUEPOS;
693 if (Angle < 5) iflag = iflag | SEFLG_TOPOCTR;
694 if (Angle == 7) Angle = 0;
695 tjd_tt = JDNDaysUT + swe_deltat_ex(JDNDaysUT, epheflag, serr);
696 Planet = DeterObject(ObjectName);
697 if (Planet != -1) {
698 if (swe_calc(tjd_tt, Planet, iflag, x, serr) == ERR)
699 return ERR;
700 } else {
701 if (call_swe_fixstar(ObjectName, tjd_tt, iflag, x, serr) == ERR)
702 return ERR;
703 }
704 if (Angle == 2 || Angle == 5) {
705 *dret = x[1];
706 } else {
707 if (Angle == 3 || Angle == 6) {
708 *dret = x[0];
709 } else {
710 xin[0] = x[0];
711 xin[1] = x[1];
712 swe_azalt(JDNDaysUT, SE_EQU2HOR, dgeo, datm[0], datm[1], xin, xaz);
713 if (Angle == 0)
714 *dret = xaz[1];
715 if (Angle == 4)
716 *dret = AppAltfromTopoAlt(xaz[1], datm[0], datm[1], helflag);
717 if (Angle == 1) {
718 xaz[0] += 180;
719 if (xaz[0] >= 360)
720 xaz[0] -= 360;
721 *dret = xaz[0];
722 }
723 }
724 }
725 return OK;
726 }
727
728 /*###################################################################
729 ' JDNDaysUT [Days]
730 ' dgeo [array: longitude, latitude, eye height above sea m]
731 ' TempE [C]
732 ' PresE [mbar]
733 ' ObjectName [-]
734 ' Angle (0 = TopoAlt, 1 = Azi, 2=Topo Declination, 3=Topo Rectascension, 4=AppAlt,5=Geo Declination, 6=Geo Rectascension)
735 ' ObjectLoc [deg]
736 */
azalt_cart(double JDNDaysUT,double * dgeo,double * datm,char * ObjectName,int32 helflag,double * dret,char * serr)737 static int32 azalt_cart(double JDNDaysUT, double *dgeo, double *datm, char *ObjectName, int32 helflag, double *dret, char *serr)
738 {
739 double x[6], xin[3], xaz[3], tjd_tt;
740 int32 Planet;
741 int32 epheflag;
742 int32 iflag = SEFLG_EQUATORIAL;
743 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
744 iflag |= epheflag;
745 if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
746 iflag |= SEFLG_NONUT | SEFLG_TRUEPOS;
747 iflag = iflag | SEFLG_TOPOCTR;
748 tjd_tt = JDNDaysUT + swe_deltat_ex(JDNDaysUT, epheflag, serr);
749 Planet = DeterObject(ObjectName);
750 if (Planet != -1) {
751 if (swe_calc(tjd_tt, Planet, iflag, x, serr) == ERR)
752 return ERR;
753 } else {
754 if (call_swe_fixstar(ObjectName, tjd_tt, iflag, x, serr) == ERR)
755 return ERR;
756 }
757 xin[0] = x[0];
758 xin[1] = x[1];
759 swe_azalt(JDNDaysUT, SE_EQU2HOR, dgeo, datm[0], datm[1], xin, xaz);
760 dret[0] = xaz[0];
761 dret[1] = xaz[1]; /* true altitude */
762 dret[2] = xaz[2]; /* apparent altitude */
763 /* also return cartesian coordinates, for apparent altitude */
764 xaz[1] = xaz[2];
765 xaz[2] = 1;
766 swi_polcart(xaz, xaz);
767 dret[3] = xaz[0];
768 dret[4] = xaz[1];
769 dret[5] = xaz[2];
770 return OK;
771 }
772
773 /*###################################################################
774 ' LatA [rad]
775 ' LongA [rad]
776 ' LatB [rad]
777 ' LongB [rad]
778 ' DistanceAngle [rad]
779 */
DistanceAngle(double LatA,double LongA,double LatB,double LongB)780 static double DistanceAngle(double LatA, double LongA, double LatB, double LongB)
781 {
782 double dlon = LongB - LongA;
783 double dlat = LatB - LatA;
784 /* Haversine formula
785 * http://www.movable-type.co.uk/scripts/GIS-FAQ-5.1.html
786 * R.W. Sinnott, Virtues of the Haversine, Sky and Telescope, vol. 68, no. 2, 1984, p. 159
787 */
788 double sindlat2 = sin(dlat / 2);
789 double sindlon2 = sin(dlon / 2);
790 double corde = sindlat2 * sindlat2 + cos(LatA) * cos(LatB) * sindlon2 *sindlon2;
791 if (corde > 1) corde = 1;
792 return 2 * asin(sqrt(corde));
793 }
794
795 /*###################################################################
796 ' heighteye [m]
797 ' TempS [C]
798 ' RH [%]
799 ' kW [-]
800 */
kW(double HeightEye,double TempS,double RH)801 static double kW(double HeightEye, double TempS, double RH)
802 {
803 /* From Schaefer , Archaeoastronomy, XV, 2000, page 128*/
804 double WT = 0.031;
805 WT *= 0.94 * (RH / 100.0) * exp(TempS / 15) * exp(-1 * HeightEye / scaleHwater);
806 return WT;
807 }
808
809 /*###################################################################
810 ' JDNDaysUT [-]
811 ' AltS [deg]
812 ' lat [deg]
813 ' kOZ [-]
814 */
kOZ(double AltS,double sunra,double Lat)815 static double kOZ(double AltS, double sunra, double Lat)
816 {
817 double CHANGEKO, OZ, LT, kOZret;
818 static TLS double koz_last, alts_last, sunra_last;
819 double altslim = 0;
820 if (AltS == alts_last && sunra == sunra_last)
821 return koz_last;
822 alts_last = AltS; sunra_last = sunra;
823 OZ = 0.031;
824 LT = Lat * DEGTORAD;
825 /* From Schaefer , Archaeoastronomy, XV, 2000, page 128*/
826 kOZret = OZ * (3.0 + 0.4 * (LT * cos(sunra * DEGTORAD) - cos(3 * LT))) / 3.0;
827 /* depending on day/night vision (altitude of sun < start astronomical twilight), KO changes from 100% to 30%
828 * see extinction section of Vistas in Astronomy page 343*/
829 altslim = -AltS - 12;
830 if (altslim < 0)
831 altslim = 0;
832 CHANGEKO = (100 - 11.6 * mymin(6, altslim)) / 100;
833 if ((0)) {
834 static int a = 0;
835 if (a == 0)
836 printf("bsk=%f %f\n", kOZret, AltS);
837 a = 1;
838 }
839 koz_last = kOZret * CHANGEKO;
840 return koz_last;
841 }
842
843 /*###################################################################
844 ' AltS [deg]
845 ' heighteye [m]
846 ' kR [-]
847 */
kR(double AltS,double HeightEye)848 static double kR(double AltS, double HeightEye)
849 {
850 /* depending on day/night vision (altitude of sun < start astronomical twilight),
851 * lambda eye sensibility changes
852 * see extinction section of Vistas in Astronomy page 343*/
853 double CHANGEK, LAMBDA;
854 double val = -AltS - 12;
855 if (val < 0) val = 0;
856 if (val > 6) val = 6;
857 /*CHANGEK = (1 - 0.166667 * Min(6, Max(-AltS - 12, 0)));*/
858 CHANGEK = (1 - 0.166667 * val );
859 LAMBDA = 0.55 + (CHANGEK - 1) * 0.04;
860 /* From Schaefer , Archaeoastronomy, XV, 2000, page 128 */
861 return 0.1066 * exp(-1 * HeightEye / scaleHrayleigh) * pow(LAMBDA / 0.55 , -4);
862 }
863
Sgn(double x)864 static int Sgn(double x)
865 {
866 if (x < 0)
867 return -1;
868 return 1;
869 }
870
871 /*###################################################################
872 ' JDNDaysUT [-]
873 ' AltS [deg]
874 ' lat [deg]
875 ' heighteye [m]
876 ' TempS [C]
877 ' RH [%]
878 ' VR [km]
879 ' ka [-]
880 */
ka(double AltS,double sunra,double Lat,double HeightEye,double TempS,double RH,double VR,char * serr)881 static double ka(double AltS, double sunra, double Lat, double HeightEye, double TempS, double RH, double VR, char *serr)
882 {
883 double CHANGEKA, LAMBDA, BetaVr, Betaa, kaact;
884 double SL = Sgn(Lat);
885 /* depending on day/night vision (altitude of sun < start astronomical twilight),
886 * lambda eye sensibility changes
887 * see extinction section of Vistas in Astronomy page 343 */
888 static TLS double alts_last, sunra_last, ka_last;
889 if (AltS == alts_last && sunra == sunra_last)
890 return ka_last;
891 alts_last = AltS; sunra_last = sunra;
892 CHANGEKA = (1 - 0.166667 * mymin(6, mymax(-AltS - 12, 0)));
893 LAMBDA = 0.55 + (CHANGEKA - 1) * 0.04;
894 if (VR != 0) {
895 if (VR >= 1) {
896 /* Visbility range from http://www1.cs.columbia.edu/CAVE/publications/pdfs/Narasimhan_CVPR03.pdf
897 * http://www.icao.int/anb/SG/AMOSSG/meetings/amossg3/wp/SN11Rev.pdf where MOR=2.995/ke
898 * factor 1.3 is the relation between "prevailing visibility" and
899 * meteorological range was derived by Koshmeider in the 1920's */
900 BetaVr = 3.912 / VR;
901 Betaa = BetaVr - (kW(HeightEye, TempS, RH) / scaleHwater + kR(AltS, HeightEye) / scaleHrayleigh) * 1000 * astr2tau;
902 kaact = Betaa * scaleHaerosol / 1000 * tau2astr;
903 if (kaact < 0) {
904 if (serr != NULL)
905 strcpy(serr, "The provided Meteorological range is too long, when taking into acount other atmospheric parameters"); /* is a warning */
906 /* return 0; * return "#HIGHVR"; */
907 }
908 } else {
909 kaact = VR - kW(HeightEye, TempS, RH) - kR(AltS, HeightEye) - kOZ(AltS, sunra, Lat);
910 if (kaact < 0) {
911 if (serr != NULL)
912 strcpy(serr, "The provided atmosphic coeefficent (ktot) is too low, when taking into acount other atmospheric parameters"); /* is a warning */
913 /* return 0; * "#LOWktot"; */
914 }
915 }
916 } else {
917 /* From Schaefer , Archaeoastronomy, XV, 2000, page 128 */
918 #ifdef SIMULATE_VICTORVB
919 if (RH <= 0.00000001) RH = 0.00000001;
920 if (RH >= 99.99999999) RH = 99.99999999;
921 #endif
922 kaact = 0.1 * exp(-1 * HeightEye / scaleHaerosol) * pow(1 - 0.32 / log(RH / 100.0), 1.33) * (1 + 0.33 * SL * sin(sunra * DEGTORAD));
923 kaact = kaact * pow(LAMBDA / 0.55, -1.3);
924 }
925 ka_last = kaact;
926 return kaact;
927 }
928
929 /*###################################################################
930 ' JDNDaysUT [-]
931 ' AltS [deg]
932 ' lat [deg]
933 ' heighteye [m]
934 ' TempS [C]
935 ' RH [%]
936 ' VR [km]
937 ' ExtType [0=ka,1=kW,2=kR,3=kOZ,4=ktot]
938 ' kt [-]
939 */
kt(double AltS,double sunra,double Lat,double HeightEye,double TempS,double RH,double VR,int32 ExtType,char * serr)940 static double kt(double AltS, double sunra, double Lat, double HeightEye, double TempS, double RH, double VR, int32 ExtType, char *serr)
941 {
942 double kRact = 0;
943 double kWact = 0;
944 double kOZact = 0;
945 double kaact = 0;
946 if (ExtType == 2 || ExtType == 4)
947 kRact = kR(AltS, HeightEye);
948 if (ExtType == 1 || ExtType == 4)
949 kWact = kW(HeightEye, TempS, RH);
950 if (ExtType == 3 || ExtType == 4)
951 kOZact = kOZ(AltS, sunra, Lat);
952 if (ExtType == 0 || ExtType == 4)
953 kaact = ka(AltS, sunra, Lat, HeightEye, TempS, RH, VR, serr);
954 if (kaact < 0)
955 kaact = 0;
956 return kWact + kRact + kOZact + kaact;
957 }
958
959 /*###################################################################
960 ' AppAlt0 [deg]
961 ' PresS [mbar]
962 ' Airmass [??]
963 */
Airmass(double AppAltO,double Press)964 static double Airmass(double AppAltO, double Press)
965 {
966 double airm, zend;
967 zend = (90 - AppAltO) * DEGTORAD;
968 if (zend > PI / 2)
969 zend = PI / 2;
970 airm = 1 / (cos(zend) + 0.025 * exp(-11 * cos(zend)));
971 return Press / 1013 * airm;
972 }
973
974 /*###################################################################
975 ' scaleH '[m]
976 ' zend [rad]
977 ' PresS [mbar]
978 ' Xext [-]
979 */
Xext(double scaleH,double zend,double Press)980 static double Xext(double scaleH, double zend, double Press)
981 {
982 return Press / 1013.0 / (cos(zend) + 0.01 * sqrt(scaleH / 1000.0) * exp(-30.0 / sqrt(scaleH / 1000.0) * cos(zend)));
983 }
984
985 /*###################################################################
986 ' scaleH '[m]
987 ' zend [rad]
988 ' PresS [mbar]
989 ' Xlay [-]
990 */
Xlay(double scaleH,double zend,double Press)991 static double Xlay(double scaleH, double zend, double Press)
992 {
993 /*return Press / 1013.0 / sqrt(1.0 - pow(sin(zend) / (1.0 + (scaleH / Ra)), 2));*/
994 double a = sin(zend) / (1.0 + (scaleH / Ra));
995 return Press / 1013.0 / sqrt(1.0 - a * a);
996 }
997
998 /*###################################################################
999 ' Meteorological formula
1000 '###################################################################
1001 ' TempS [C]
1002 ' HeightEye [m]
1003 ' TempEfromTempS [C]
1004 */
TempEfromTempS(double TempS,double HeightEye,double Lapse)1005 static double TempEfromTempS(double TempS, double HeightEye, double Lapse)
1006 {
1007 return TempS - Lapse * HeightEye;
1008 }
1009
1010 /*###################################################################
1011 ' TempS [C]
1012 ' PresS [mbar]
1013 ' HeightEye [m]
1014 ' PresEfromPresS [mbar]
1015 */
PresEfromPresS(double TempS,double Press,double HeightEye)1016 static double PresEfromPresS(double TempS, double Press, double HeightEye)
1017 {
1018 return Press * exp(-9.80665 * 0.0289644 / (Kelvin(TempS) + 3.25 * HeightEye / 1000) / 8.31441 * HeightEye);
1019 }
1020
1021 /*###################################################################
1022 ' AltO [deg]
1023 ' JDNDaysUT [-]
1024 ' AltS [deg]
1025 ' lat [deg]
1026 ' heighteye [m]
1027 ' TempS [C]
1028 ' PresS [mbar]
1029 ' RH [%]
1030 ' VR [km]
1031 ' Deltam [-]
1032 */
Deltam(double AltO,double AltS,double sunra,double Lat,double HeightEye,double * datm,int32 helflag,char * serr)1033 static double Deltam(double AltO, double AltS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, char *serr)
1034 {
1035 double zend, xR, XW, Xa, XOZ;
1036 double PresE = PresEfromPresS(datm[1], datm[0], HeightEye);
1037 double TempE = TempEfromTempS(datm[1], HeightEye, LapseSA);
1038 double AppAltO = AppAltfromTopoAlt(AltO, TempE, PresE, helflag);
1039 double deltam;
1040 static TLS double alts_last, alto_last, sunra_last, deltam_last;
1041 if (AltS == alts_last && AltO == alto_last && sunra == sunra_last)
1042 return deltam_last;
1043 alts_last = AltS; alto_last = AltO; sunra_last = sunra;
1044 if (staticAirmass == 0) {
1045 zend = (90 - AppAltO) * DEGTORAD;
1046 if (zend > PI / 2)
1047 zend = PI / 2;
1048 /* From Schaefer , Archaeoastronomy, XV, 2000, page 128*/
1049 xR = Xext(scaleHrayleigh, zend, datm[0]);
1050 XW = Xext(scaleHwater, zend, datm[0]);
1051 Xa = Xext(scaleHaerosol, zend, datm[0]);
1052 XOZ = Xlay(scaleHozone, zend, datm[0]);
1053 deltam = kR(AltS, HeightEye) * xR + kt(AltS, sunra, Lat, HeightEye, datm[1], datm[2], datm[3], 0, serr) * Xa + kOZ(AltS, sunra, Lat) * XOZ + kW(HeightEye, datm[1], datm[2]) * XW;
1054 } else {
1055 deltam = kt(AltS, sunra, Lat, HeightEye, datm[1], datm[2], datm[3], 4, serr) * Airmass(AppAltO, datm[0]);
1056 }
1057 deltam_last = deltam;
1058 return deltam;
1059 }
1060
1061 /*###################################################################
1062 ' AltO [deg]
1063 ' JDNDaysUT [-]
1064 ' AltS [deg]
1065 ' lat [deg]
1066 ' heighteye [m]
1067 ' TempS [C]
1068 ' PresS [mbar]
1069 ' RH [%]
1070 ' VR [km]
1071 ' Bn [nL]
1072 */
Bn(double AltO,double JDNDayUT,double AltS,double sunra,double Lat,double HeightEye,double * datm,int32 helflag,char * serr)1073 static double Bn(double AltO, double JDNDayUT, double AltS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, char *serr)
1074 {
1075 double PresE = PresEfromPresS(datm[1], datm[0], HeightEye);
1076 double TempE = TempEfromTempS(datm[1], HeightEye, LapseSA);
1077 double AppAltO = AppAltfromTopoAlt(AltO, TempE, PresE, helflag);
1078 double zend, YearB, MonthB, DayB, Bna, kX, Bnb;
1079 double B0 = 0.0000000000001, dut;
1080 int iyar, imon, iday;
1081 /* Below altitude of 10 degrees, the Bn stays the same (see page 343 Vistas in Astronomy) */
1082 if (AppAltO < 10)
1083 AppAltO = 10;
1084 zend = (90 - AppAltO) * DEGTORAD;
1085 /* From Schaefer , Archaeoastronomy, XV, 2000, page 128 and adjusted for sunspot period*/
1086 /*YearB = DatefromJDut(JDNDayUT, 1);
1087 MonthB = DatefromJDut(JDNDayUT, 2);
1088 DayB = DatefromJDut(JDNDayUT, 3);*/
1089 swe_revjul(JDNDayUT, SE_GREG_CAL, &iyar, &imon, &iday, &dut);
1090 YearB = iyar; MonthB = imon; DayB = iday;
1091 Bna = B0 * (1 + 0.3 * cos(6.283 * (YearB + ((DayB - 1) / 30.4 + MonthB - 1) / 12 - 1990.33) / 11.1));
1092 kX = Deltam(AltO, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
1093 /* From Schaefer , Archaeoastronomy, XV, 2000, page 129 */
1094 Bnb = Bna * (0.4 + 0.6 / sqrt(1 - 0.96 * pow(sin(zend), 2))) * pow(10, -0.4 * kX);
1095 return mymax(Bnb, 0) * erg2nL;
1096 }
1097
1098 /*###################################################################
1099 ' JDNDaysUT [-]
1100 ' dgeo [array: longitude, latitude, eye height above sea m]
1101 ' TempE [C]
1102 ' PresE [mbar]
1103 ' ObjectName [-]
1104 ' Magnitude [-]
1105 */
Magnitude(double JDNDaysUT,double * dgeo,char * ObjectName,int32 helflag,double * dmag,char * serr)1106 static int32 Magnitude(double JDNDaysUT, double *dgeo, char *ObjectName, int32 helflag, double *dmag, char *serr)
1107 {
1108 double x[20];
1109 int32 Planet, iflag, epheflag;
1110 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
1111 *dmag = -99.0;
1112 Planet = DeterObject(ObjectName);
1113 iflag = SEFLG_TOPOCTR | SEFLG_EQUATORIAL | epheflag;
1114 if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
1115 iflag |= SEFLG_NONUT|SEFLG_TRUEPOS;
1116 if (Planet != -1) {
1117 /**dmag = Phenomena(JDNDaysUT, Lat, Longitude, HeightEye, TempE, PresE, ObjectName, 4);*/
1118 swe_set_topo(dgeo[0], dgeo[1], dgeo[2]);
1119 if (swe_pheno_ut(JDNDaysUT, Planet, iflag, x, serr) == ERR)
1120 return ERR;
1121 *dmag = x[4];
1122 } else {
1123 if (call_swe_fixstar_mag(ObjectName, dmag, serr) == ERR)
1124 return ERR;
1125 }
1126 return OK;
1127 }
1128
1129 #if 0
1130 static int32 fast_magnitude(double tjd, double *dgeo, char *ObjectName, int32 helflag, double *dmag, char *serr)
1131 {
1132 int32 retval = OK, ipl, ipli;
1133 double dtjd;
1134 static TLS double tjdsv[3];
1135 static TLS double dmagsv[3];
1136 static TLS int32 helflagsv[3];
1137 ipl = DeterObject(ObjectName);
1138 ipli = ipl;
1139 if (ipli > SE_MOON)
1140 ipli = 2;
1141 dtjd = tjd - tjdsv[ipli];
1142 if (tjdsv[ipli] != 0 && helflag == helflagsv[ipli] && fabs(dtjd) < 5.0 / 1440.0) {
1143 *dmag = dmagsv[ipli];
1144 } else {
1145 retval = Magnitude(tjd, dgeo, ObjectName, helflag, dmag, serr);
1146 tjdsv[ipli] = tjd;
1147 helflagsv[ipli] = helflag;
1148 dmagsv[ipli] = *dmag;
1149 }
1150 return retval;
1151 }
1152 #endif
1153
1154 /*###################################################################
1155 ' dist [km]
1156 ' phasemoon [-]
1157 ' MoonsBrightness [-]
1158 */
MoonsBrightness(double dist,double phasemoon)1159 static double MoonsBrightness(double dist, double phasemoon)
1160 {
1161 double log10 = 2.302585092994;
1162 /*Moon's brightness changes with distance: http://hem.passagen.se/pausch/comp/ppcomp.html#15 */
1163 return -21.62 + 5 * log(dist / (Ra / 1000)) / log10 + 0.026 * fabs(phasemoon) + 0.000000004 * pow(phasemoon, 4);
1164 }
1165
1166 /*###################################################################
1167 ' AltM [deg]
1168 ' AziM [deg]
1169 ' AziS [deg]
1170 ' MoonPhase [deg]
1171 */
MoonPhase(double AltM,double AziM,double AltS,double AziS)1172 static double MoonPhase(double AltM, double AziM, double AltS, double AziS)
1173 {
1174 double AltMi = AltM * DEGTORAD;
1175 double AltSi = AltS * DEGTORAD;
1176 double AziMi = AziM * DEGTORAD;
1177 double AziSi = AziS * DEGTORAD;
1178 double MoonAvgPar = 0.95;
1179 // return 180 - acos(cos(AziSi - AziMi) * cos(AltMi + MoonAvgPar * DEGTORAD) * cos(AltSi) + sin(AltSi) * sin(AltMi + MoonAvgPar * DEGTORAD)) / DEGTORAD;
1180 return 180 - acos(cos(AziSi - AziMi - MoonAvgPar * DEGTORAD) * cos(AltMi + MoonAvgPar * DEGTORAD) * cos(AltSi) + sin(AltSi) * sin(AltMi + MoonAvgPar * DEGTORAD)) / DEGTORAD;
1181 }
1182
1183 /*###################################################################
1184 ' Pressure [mbar]
1185 */
Bm(double AltO,double AziO,double AltM,double AziM,double AltS,double AziS,double sunra,double Lat,double HeightEye,double * datm,int32 helflag,char * serr)1186 static double Bm(double AltO, double AziO, double AltM, double AziM, double AltS, double AziS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, char *serr)
1187 {
1188 double M0 = -11.05;
1189 double Bm = 0;
1190 double RM, kXM, kX, C3, FM, phasemoon, MM;
1191 double lunar_radius = 0.25 * DEGTORAD;
1192 AS_BOOL object_is_moon = FALSE;
1193 if (AltO == AltM && AziO == AziM)
1194 object_is_moon = TRUE;
1195 if (AltM > -0.26 && !object_is_moon) { // second condition added by Dieter, SE2.06
1196 /* moon only adds light when (partly) above horizon
1197 * From Schaefer , Archaeoastronomy, XV, 2000, page 129*/
1198 RM = DistanceAngle(AltO * DEGTORAD, AziO * DEGTORAD, AltM * DEGTORAD, AziM * DEGTORAD) / DEGTORAD;
1199 if (RM <= lunar_radius) // addition by Dieter for objects behind the Moon, SE2.06
1200 RM = lunar_radius;
1201 kXM = Deltam(AltM, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
1202 kX = Deltam(AltO, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
1203 C3 = pow(10, -0.4 * kXM);
1204 FM = (62000000.0) / RM / RM + pow(10, 6.15 - RM / 40) + pow(10, 5.36) * (1.06 + pow(cos(RM * DEGTORAD), 2));
1205 Bm = FM * C3 + 440000 * (1 - C3);
1206 phasemoon = MoonPhase(AltM, AziM, AltS, AziS);
1207 MM = MoonsBrightness(MoonDistance, phasemoon);
1208 Bm = Bm * pow(10, -0.4 * (MM - M0 + 43.27));
1209 Bm = Bm * (1 - pow(10, -0.4 * kX));
1210 }
1211 Bm = mymax(Bm, 0) * erg2nL;
1212 return Bm;
1213 }
1214
1215 /*###################################################################
1216 ' Pressure [mbar]
1217 */
Btwi(double AltO,double AziO,double AltS,double AziS,double sunra,double Lat,double HeightEye,double * datm,int32 helflag,char * serr)1218 static double Btwi(double AltO, double AziO, double AltS, double AziS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, char *serr)
1219 {
1220 double M0 = -11.05;
1221 double MS = -26.74;
1222 double PresE = PresEfromPresS(datm[1], datm[0], HeightEye);
1223 double TempE = TempEfromTempS(datm[1], HeightEye, LapseSA);
1224 double AppAltO = AppAltfromTopoAlt(AltO, TempE, PresE, helflag);
1225 double ZendO = 90 - AppAltO;
1226 double RS = DistanceAngle(AltO * DEGTORAD, AziO * DEGTORAD, AltS * DEGTORAD, AziS * DEGTORAD) / DEGTORAD;
1227 double kX = Deltam(AltO, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
1228 double k = kt(AltS, sunra, Lat, HeightEye, datm[1], datm[2], datm[3], 4, serr);
1229 /* From Schaefer , Archaeoastronomy, XV, 2000, page 129*/
1230 double Btwi = pow(10, -0.4 * (MS - M0 + 32.5 - AltS - (ZendO / (360 * k))));
1231 Btwi = Btwi * (100 / RS) * (1 - pow(10, -0.4 * kX));
1232 Btwi = mymax(Btwi, 0) * erg2nL;
1233 return Btwi;
1234 }
1235
1236 /*###################################################################
1237 ' Pressure [mbar]
1238 2300 REM Daylight brightness
1239 2310 C4=10.0^(-.4*K(I)*XS)
1240 2320 FS=6.2E+07*(RS^-2)+(10^(6.15-RS/40))
1241 2330 FS=FS+(10^5.36)*(1.06+((COS(RS*RD))^2))
1242 2340 BD=10^(-.4*(MS(I)-MO(I)+43.27))
1243 2350 BD=BD*(1-10^(-.4*K(I)*X))
1244 2360 BD=BD*(FS*C4+440000.0*(1-C4))
1245 */
Bday(double AltO,double AziO,double AltS,double AziS,double sunra,double Lat,double HeightEye,double * datm,int32 helflag,char * serr)1246 static double Bday(double AltO, double AziO, double AltS, double AziS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, char *serr)
1247 {
1248 double M0 = -11.05;
1249 double MS = -26.74;
1250 double RS = DistanceAngle(AltO * DEGTORAD, AziO * DEGTORAD, AltS * DEGTORAD, AziS * DEGTORAD) / DEGTORAD;
1251 double kXS = Deltam(AltS, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
1252 double kX = Deltam(AltO, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
1253 /* From Schaefer , Archaeoastronomy, XV, 2000, page 129*/
1254 double C4 = pow(10, -0.4 * kXS);
1255 double FS = (62000000.0) / RS / RS + pow(10, (6.15 - RS / 40)) + pow(10, 5.36) * (1.06 + pow(cos(RS * DEGTORAD), 2));
1256 double Bday = FS * C4 + 440000.0 * (1 - C4);
1257 Bday = Bday * pow(10, (-0.4 * (MS - M0 + 43.27)));
1258 Bday = Bday * (1 - pow(10, -0.4 * kX));
1259 Bday = mymax(Bday, 0) * erg2nL;
1260 return Bday;
1261 }
1262
1263 /*###################################################################
1264 ' Value [nL]
1265 ' PresS [mbar]
1266 ' Bcity [nL]
1267 */
Bcity(double Value,double Press)1268 static double Bcity(double Value, double Press)
1269 {
1270 double Bcity = Value;
1271 Press += 0.0; /* unused; statement prevents compiler warning */
1272 Bcity = mymax(Bcity, 0);
1273 return Bcity;
1274 }
1275
1276 /*###################################################################
1277 ' Pressure [mbar]
1278 */
Bsky(double AltO,double AziO,double AltM,double AziM,double JDNDaysUT,double AltS,double AziS,double sunra,double Lat,double HeightEye,double * datm,int32 helflag,char * serr)1279 static double Bsky(double AltO, double AziO, double AltM, double AziM, double JDNDaysUT, double AltS, double AziS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, char *serr)
1280 {
1281 double Bsky = 0;
1282 if (AltS < -3) {
1283 Bsky += Btwi(AltO, AziO, AltS, AziS, sunra, Lat, HeightEye, datm, helflag, serr);
1284 } else {
1285 if (AltS > 4) {
1286 Bsky += Bday(AltO, AziO, AltS, AziS, sunra, Lat, HeightEye, datm, helflag, serr);
1287 } else {
1288 Bsky += mymin(Bday(AltO, AziO, AltS, AziS, sunra, Lat, HeightEye, datm, helflag, serr), Btwi(AltO, AziO, AltS, AziS, sunra, Lat, HeightEye, datm, helflag, serr));
1289 if ((0)) {
1290 static int a = 0;
1291 if (a == 0)
1292 printf("bsk=%f\n", Bsky);
1293 a = 1;
1294 }
1295 }
1296 }
1297 /* if max. Bm [1E7] <5% of Bsky don't add Bm*/
1298 if (Bsky < 200000000.0)
1299 Bsky += Bm(AltO, AziO, AltM, AziM, AltS, AziS, sunra, Lat, HeightEye, datm, helflag, serr);
1300 if (AltS <= 0)
1301 Bsky += Bcity(0, datm[0]);
1302 /* if max. Bn [250] <5% of Bsky don't add Bn*/
1303 if (Bsky < 5000)
1304 Bsky = Bsky + Bn(AltO, JDNDaysUT, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
1305 /* if max. Bm [1E7] <5% of Bsky don't add Bm*/
1306 return Bsky;
1307 }
1308
1309 /* default handling:
1310 * 1. datm (atmospheric conditions):
1311 * datm consists of
1312 * [0] atmospheric pressure
1313 * [1] temperature
1314 * [2] relative humidity
1315 * [3] extinction coefficient
1316 * In order to get default values for [0..2], set datm[0] = 0.
1317 * Default values for [1-2] are only provided if [0] == 0.
1318 * [3] defaults outside this function, depending on [0-2].
1319 *
1320 * 2. dobs (observer definition):
1321 * [0] age (default 36)
1322 * [1] Snellen ratio or visual acuity of observer (default 1)
1323 */
default_heliacal_parameters(double * datm,double * dgeo,double * dobs,int helflag)1324 static void default_heliacal_parameters(double *datm, double *dgeo, double *dobs, int helflag)
1325 {
1326 int i;
1327 if (datm[0] <= 0) {
1328 /* estimate atmospheric pressure, according to the
1329 * International Standard Atmosphere (ISA) */
1330 datm[0] = 1013.25 * pow(1 - 0.0065 * dgeo[2] / 288, 5.255);
1331 /* temperature */
1332 if (datm[1] == 0)
1333 datm[1] = 15 - 0.0065 * dgeo[2];
1334 /* relative humidity, independent of atmospheric pressure and altitude */
1335 if (datm[2] == 0)
1336 datm[2] = 40;
1337 /* note: datm[3] / VR defaults outside this function */
1338 } else {
1339 #ifndef SIMULATE_VICTORVB
1340 if (datm[2] <= 0.00000001) datm[2] = 0.00000001;
1341 if (datm[2] >= 99.99999999) datm[2] = 99.99999999;
1342 #endif
1343 }
1344 /* age of observer */
1345 if (dobs[0] == 0)
1346 dobs[0] = 36;
1347 /* SN Snellen factor of the visual acuity of the observer */
1348 if (dobs[1] == 0)
1349 dobs[1] = 1;
1350 if (!(helflag & SE_HELFLAG_OPTICAL_PARAMS)) {
1351 for (i = 2; i <= 5; i++)
1352 dobs[i] = 0;
1353 }
1354 /* OpticMagn undefined -> use eye */
1355 if (dobs[3] == 0) {
1356 dobs[2] = 1; /* Binocular = 1 */
1357 dobs[3] = 1; /* OpticMagn = 1: use eye */
1358 /* dobs[4] and dobs[5] (OpticDia and OpticTrans) will be defaulted in
1359 * OpticFactor() */
1360 }
1361 }
1362
1363 /*###################################################################
1364 ' age [Year]
1365 ' SN [-]
1366 ' AltO [deg]
1367 ' AziO [deg]
1368 ' AltM [deg]
1369 ' AziM [deg]
1370 ' MoonDistance [km]
1371 ' JDNDaysUT [-]
1372 ' AltS [deg]
1373 ' AziS [deg]
1374 ' lat [deg]
1375 ' heighteye [m]
1376 ' TempS [C]
1377 ' PresS [mbar]
1378 ' RH [%]
1379 ' VR [km]
1380 ' VisLimMagn [-]
1381 */
VisLimMagn(double * dobs,double AltO,double AziO,double AltM,double AziM,double JDNDaysUT,double AltS,double AziS,double sunra,double Lat,double HeightEye,double * datm,int32 helflag,int32 * scotopic_flag,char * serr)1382 static double VisLimMagn(double *dobs, double AltO, double AziO, double AltM, double AziM, double JDNDaysUT, double AltS, double AziS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, int32 *scotopic_flag, char *serr)
1383 {
1384 double C1, C2, Th, kX, Bsk, CorrFactor1, CorrFactor2;
1385 double log10 = 2.302585092994;
1386 AS_BOOL is_scotopic = FALSE;
1387 /*double Age = dobs[0];*/
1388 /*double SN = dobs[1];*/
1389 Bsk = Bsky(AltO, AziO, AltM, AziM, JDNDaysUT, AltS, AziS, sunra, Lat, HeightEye, datm, helflag, serr);
1390 /* Schaefer, Astronomy and the limits of vision, Archaeoastronomy, 1993 Verder:*/
1391 kX = Deltam(AltO, AltS, sunra, Lat, HeightEye, datm, helflag, serr);
1392 if ((0)) {
1393 static int a = 0;
1394 if (a == 0)
1395 printf("bsk=%f, kx=%f\n", Bsk, kX);
1396 a = 1;
1397 }
1398 /* influence of age*/
1399 /*Fa = mymax(1, pow(p(23, Bsk) / p(Age, Bsk), 2)); */
1400 CorrFactor1 = OpticFactor(Bsk, kX, dobs, JDNDaysUT, "", 1, helflag);
1401 CorrFactor2 = OpticFactor(Bsk, kX, dobs, JDNDaysUT, "", 0, helflag);
1402 //if (Bsk < BNIGHT)
1403 if (Bsk < 1645) /* use this function for BNIGHT to make the function continuous */
1404 is_scotopic = TRUE;
1405 if (helflag & SE_HELFLAG_VISLIM_PHOTOPIC)
1406 is_scotopic = FALSE;
1407 if (helflag & SE_HELFLAG_VISLIM_SCOTOPIC)
1408 is_scotopic = TRUE;
1409 /* From Schaefer , Archaeoastronomy, XV, 2000, page 129*/
1410 if (is_scotopic) {
1411 C1 = 1.5848931924611e-10; /*pow(10, -9.8);*/ /* C1 = 10 ^ (-9.8);*/
1412 C2 = 0.012589254117942; /*pow(10, -1.9);*/ /* C2 = 10 ^ (-1.9);*/
1413 if (scotopic_flag != NULL)
1414 *scotopic_flag = 1;
1415 } else {
1416 C1 = 4.4668359215096e-9; /*pow(10, -8.35);*/ /* C1 = 10 ^ (-8.35);*/
1417 C2 = 1.2589254117942e-6; /*pow(10, -5.9);*/ /* C2 = 10 ^ (-5.9);*/
1418 if (scotopic_flag != NULL)
1419 *scotopic_flag = 0;
1420 }
1421 if (scotopic_flag != NULL) {
1422 if (BNIGHT * BNIGHT_FACTOR > Bsk && BNIGHT / BNIGHT_FACTOR < Bsk)
1423 *scotopic_flag |= 2;
1424 }
1425 /*Th = C1 * pow(1 + sqrt(C2 * Bsk), 2) * Fa;*/
1426 /*Bsk = Bsk / CorrFactor1;*/
1427 Bsk = Bsk * CorrFactor1;
1428 Th = C1 * pow(1 + sqrt(C2 * Bsk), 2) * CorrFactor2;
1429 #if SWEHEL_DEBUG
1430 fprintf(stderr, "Bsk=%f, ", Bsk);
1431 fprintf(stderr, "kX =%f, ", kX);
1432 fprintf(stderr, "Th =%f, ", Th);
1433 fprintf(stderr, "CorrFactor1=%f, ", CorrFactor1);
1434 fprintf(stderr, "CorrFactor2=%f\n", CorrFactor2);
1435 #endif
1436 /* Visual limiting magnitude of point source*/
1437 #if 0
1438 if (SN <= 0.00000001)
1439 SN = 0.00000001;
1440 return -16.57 - 2.5 * (log(Th) / log10) - kX + 5.0 * (log(SN) / log10);*/
1441 #endif
1442 return -16.57 - 2.5 * (log(Th) / log10);
1443 }
1444
1445 /* tolower star name, but not Bayer designation */
tolower_string_star(char * str)1446 static char *tolower_string_star(char *str)
1447 {
1448 char *sp;
1449 for (sp = str; *sp != '\0' && *sp != ','; sp++)
1450 *sp = tolower(*sp);
1451 return str;
1452 }
1453
1454 /* Limiting magnitude in dark skies
1455 * for information about input parameters, see function swe_heliacal_ut().
1456 *
1457 * function returns:
1458 * -1 Error
1459 * -2 Object is below horizon
1460 * 0 OK, photopic vision
1461 * |1 OK, scotopic vision
1462 * |2 OK, near limit photopic/scotopic
1463 */
swe_vis_limit_mag(double tjdut,double * dgeo,double * datm,double * dobs,char * ObjectName,int32 helflag,double * dret,char * serr)1464 int32 CALL_CONV swe_vis_limit_mag(double tjdut, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 helflag, double *dret, char *serr)
1465 {
1466 int32 retval = OK, i, scotopic_flag = 0;
1467 double AltO, AziO, AltM, AziM, AltS, AziS;
1468 double sunra;
1469 for (i = 0; i < 7; i++)
1470 dret[i] = 0;
1471 tolower_string_star(ObjectName);
1472 if (DeterObject(ObjectName) == SE_SUN) {
1473 if (serr != NULL) {
1474 strcpy(serr, "it makes no sense to call swe_vis_limit_mag() for the Sun");
1475 }
1476 return ERR;
1477 }
1478 swi_set_tid_acc(tjdut, helflag, 0, serr);
1479 sunra = SunRA(tjdut, helflag, serr);
1480 default_heliacal_parameters(datm, dgeo, dobs, helflag);
1481 swe_set_topo(dgeo[0], dgeo[1], dgeo[2]);
1482 if (ObjectLoc(tjdut, dgeo, datm, ObjectName, 0, helflag, &AltO, serr) == ERR)
1483 return ERR;
1484 if (AltO < 0) {
1485 if (serr != NULL)
1486 strcpy(serr, "object is below local horizon");
1487 *dret = -100;
1488 return -2;
1489 }
1490 if (ObjectLoc(tjdut, dgeo, datm, ObjectName, 1, helflag, &AziO, serr) == ERR)
1491 return ERR;
1492 if (helflag & SE_HELFLAG_VISLIM_DARK) {
1493 AltS = -90;
1494 AziS = 0;
1495 } else {
1496 if (ObjectLoc(tjdut, dgeo, datm, "sun", 0, helflag, &AltS, serr) == ERR)
1497 return ERR;
1498 if (ObjectLoc(tjdut, dgeo, datm, "sun", 1, helflag, &AziS, serr) == ERR)
1499 return ERR;
1500 }
1501 if (strncmp(ObjectName, "moon", 4) == 0 ||
1502 (helflag & SE_HELFLAG_VISLIM_DARK) ||
1503 (helflag & SE_HELFLAG_VISLIM_NOMOON)
1504 ) {
1505 AltM = -90; AziM = 0;
1506 } else {
1507 if (ObjectLoc(tjdut, dgeo, datm, "moon", 0, helflag, &AltM, serr) == ERR)
1508 return ERR;
1509 if (ObjectLoc(tjdut, dgeo, datm, "moon", 1, helflag, &AziM, serr) == ERR)
1510 return ERR;
1511 }
1512 #if SWEHEL_DEBUG
1513 {
1514 int i;
1515 for (i = 0; i < 6;i++)
1516 printf("dobs[%d] = %f\n", i, dobs[i]);
1517 printf("AltO = %.10f, AziO = %.10f\n", AltO, AziO);
1518 printf("AltM = %.10f, AziM = %.10f\n", AltM, AziM);
1519 printf("AltS = %.10f, AziS = %.10f\n", AltS, AziS);
1520 printf("JD = %.10f\n", tjdut);
1521 printf("lat = %f, eyeh = %f\n", dgeo[1], dgeo[2]);
1522 for (i = 0; i < 4;i++)
1523 printf("datm[%d] = %f\n", i, datm[i]);
1524 printf("helflag = %d\n", helflag);
1525 }
1526 #endif
1527 dret[0] = VisLimMagn(dobs, AltO, AziO, AltM, AziM, tjdut, AltS, AziS, sunra, dgeo[1], dgeo[2], datm, helflag, &scotopic_flag, serr);
1528 dret[1] = AltO;
1529 dret[2] = AziO;
1530 dret[3] = AltS;
1531 dret[4] = AziS;
1532 dret[5] = AltM;
1533 dret[6] = AziM;
1534 if (Magnitude(tjdut, dgeo, ObjectName, helflag, &(dret[7]), serr) == ERR)
1535 return ERR;
1536 retval = scotopic_flag;
1537 /*dret[8] = (double) is_scotopic;*/
1538 /*if (*serr != '\0') * in VisLimMagn(), serr is only a warning *
1539 retval = ERR; */
1540 return retval;
1541 }
1542
1543 /*###################################################################
1544 ' Magn [-]
1545 ' age [Year]
1546 ' SN [-]
1547 ' AltO [deg]
1548 ' AziO [deg]
1549 ' AltM [deg]
1550 ' AziM [deg]
1551 ' MoonDistance [km]
1552 ' JDNDaysUT [-]
1553 ' AziS [deg]
1554 ' lat [deg]
1555 ' heighteye [m]
1556 ' Temperature [C]
1557 ' Pressure [mbar]
1558 ' RH [%]
1559 ' VR [km]
1560 ' TopoArcVisionis [deg]
1561 */
TopoArcVisionis(double Magn,double * dobs,double AltO,double AziO,double AltM,double AziM,double JDNDaysUT,double AziS,double sunra,double Lat,double HeightEye,double * datm,int32 helflag,double * dret,char * serr)1562 static int32 TopoArcVisionis(double Magn, double *dobs, double AltO, double AziO, double AltM, double AziM, double JDNDaysUT, double AziS, double sunra, double Lat, double HeightEye, double *datm, int32 helflag, double *dret, char *serr)
1563 {
1564 double Xm, Ym, AltSi, AziSi;
1565 double xR = 0;
1566 double Xl = 45;
1567 double Yl, Yr;
1568 Yl = Magn - VisLimMagn(dobs, AltO, AziO, AltM, AziM, JDNDaysUT, AltO - Xl, AziS, sunra, Lat, HeightEye, datm, helflag, NULL, serr);
1569 /* if (*serr != '\0') return ERR; * serr is only a warning */
1570 Yr = Magn - VisLimMagn(dobs, AltO, AziO, AltM, AziM, JDNDaysUT, AltO - xR, AziS, sunra, Lat, HeightEye, datm, helflag, NULL, serr);
1571 /* if (*serr != '\0') return ERR; * serr is only a warning */
1572 /* http://en.wikipedia.org/wiki/Bisection_method*/
1573 if ((Yl * Yr) <= 0) {
1574 while(fabs(xR - Xl) > epsilon) {
1575 /*Calculate midpoint of domain*/
1576 Xm = (xR + Xl) / 2.0;
1577 AltSi = AltO - Xm;
1578 AziSi = AziS;
1579 Ym = Magn - VisLimMagn(dobs, AltO, AziO, AltM, AziM, JDNDaysUT, AltSi, AziSi, sunra, Lat, HeightEye, datm, helflag, NULL, serr);
1580 /* if (*serr != '\0') return ERR; * serr is only a warning */
1581 if ((Yl * Ym) > 0) {
1582 /* Throw away left half*/
1583 Xl = Xm;
1584 Yl = Ym;
1585 } else {
1586 /* Throw away right half */
1587 xR = Xm;
1588 Yr = Ym;
1589 }
1590 }
1591 Xm = (xR + Xl) / 2.0;
1592 } else {
1593 Xm = 99;
1594 }
1595 if (Xm < AltO)
1596 Xm = AltO;
1597 *dret = Xm;
1598 return OK;
1599 }
1600
swe_topo_arcus_visionis(double tjdut,double * dgeo,double * datm,double * dobs,int32 helflag,double mag,double azi_obj,double alt_obj,double azi_sun,double azi_moon,double alt_moon,double * dret,char * serr)1601 int32 CALL_CONV swe_topo_arcus_visionis(double tjdut, double *dgeo, double *datm, double *dobs, int32 helflag, double mag, double azi_obj, double alt_obj, double azi_sun, double azi_moon, double alt_moon, double *dret, char *serr)
1602 {
1603 double sunra;
1604 swi_set_tid_acc(tjdut, helflag, 0, serr);
1605 sunra = SunRA(tjdut, helflag, serr);
1606 if (serr != NULL && *serr != '\0')
1607 return ERR;
1608 default_heliacal_parameters(datm, dgeo, dobs, helflag);
1609 return TopoArcVisionis(mag, dobs, alt_obj, azi_obj, alt_moon, azi_moon, tjdut, azi_sun, sunra, dgeo[1], dgeo[2], datm, helflag, dret, serr);
1610 }
1611
1612 /*###################################################################*/
1613 /*' Magn [-]
1614 ' age [Year]
1615 ' SN Snellen factor of the visual aquity of the observer
1616 see: http://www.i-see.org/eyecharts.html#make-your-own
1617 ' AziO [deg]
1618 ' AltM [deg]
1619 ' AziM [deg]
1620 ' MoonDistance [km]
1621 ' JDNDaysUT [-]
1622 ' AziS [deg]
1623 ' Lat [deg]
1624 ' HeightEye [m]
1625 ' Temperature [C]
1626 ' Pressure [mbar]
1627 ' RH [%] relative humidity
1628 ' VR [km] Meteorological Range,
1629 see http://www.iol.ie/~geniet/eng/atmoastroextinction.htm
1630 ' TypeAngle
1631 ' [0=Object's altitude,
1632 ' 1=Arcus Visonis (Object's altitude - Sun's altitude),
1633 ' 2=Sun's altitude]
1634 ' HeliacalAngle [deg]
1635 */
HeliacalAngle(double Magn,double * dobs,double AziO,double AltM,double AziM,double JDNDaysUT,double AziS,double * dgeo,double * datm,int32 helflag,double * dangret,char * serr)1636 static int32 HeliacalAngle(double Magn, double *dobs, double AziO, double AltM, double AziM, double JDNDaysUT, double AziS, double *dgeo, double *datm, int32 helflag, double *dangret, char *serr)
1637 {
1638 double x, minx, maxx, xmin, ymin, Xl, xR, Yr, Yl, Xm, Ym, xmd, ymd;
1639 double Arc, DELTAx;
1640 double sunra = SunRA(JDNDaysUT, helflag, serr);
1641 double Lat = dgeo[1];
1642 double HeightEye = dgeo[2];
1643 if (PLSV == 1) {
1644 dangret[0] = criticalangle;
1645 dangret[1] = criticalangle + Magn * 2.492 + 13.447;
1646 dangret[2] = -(Magn * 2.492 + 13.447); /* Magn * 1.1 + 8.9;*/
1647 return OK;
1648 }
1649 minx = 2;
1650 maxx = 20;
1651 xmin = 0;
1652 ymin = 10000;
1653 for (x = minx; x <= maxx; x++) {
1654 if (TopoArcVisionis(Magn, dobs, x, AziO, AltM, AziM, JDNDaysUT, AziS, sunra, Lat, HeightEye, datm, helflag, &Arc, serr) == ERR)
1655 return ERR;
1656 if (Arc < ymin) {
1657 ymin = Arc;
1658 xmin = x;
1659 }
1660 }
1661 Xl = xmin - 1;
1662 xR = xmin + 1;
1663 if (TopoArcVisionis(Magn, dobs, xR, AziO, AltM, AziM, JDNDaysUT, AziS, sunra, Lat, HeightEye, datm, helflag, &Yr, serr) == ERR)
1664 return ERR;
1665 if (TopoArcVisionis(Magn, dobs, Xl, AziO, AltM, AziM, JDNDaysUT, AziS, sunra, Lat, HeightEye, datm, helflag, &Yl, serr) == ERR)
1666 return ERR;
1667 /* http://en.wikipedia.org/wiki/Bisection_method*/
1668 while(fabs(xR - Xl) > 0.1) {
1669 /* Calculate midpoint of domain */
1670 Xm = (xR + Xl) / 2.0;
1671 DELTAx = 0.025;
1672 xmd = Xm + DELTAx;
1673 if (TopoArcVisionis(Magn, dobs, Xm, AziO, AltM, AziM, JDNDaysUT, AziS, sunra, Lat, HeightEye, datm, helflag, &Ym, serr) == ERR)
1674 return ERR;
1675 if (TopoArcVisionis(Magn, dobs, xmd, AziO, AltM, AziM, JDNDaysUT, AziS, sunra, Lat, HeightEye, datm, helflag, &ymd, serr) == ERR)
1676 return ERR;
1677 if (Ym >= ymd) {
1678 /* Throw away left half */
1679 Xl = Xm;
1680 Yl = Ym;
1681 } else {
1682 /*Throw away right half */
1683 xR = Xm;
1684 Yr = Ym;
1685 }
1686 }
1687 Xm = (xR + Xl) / 2.0;
1688 Ym = (Yr + Yl) / 2.0;
1689 dangret[1] = Ym;
1690 dangret[2] = Xm - Ym;
1691 dangret[0] = Xm;
1692 return OK;
1693 }
1694
swe_heliacal_angle(double tjdut,double * dgeo,double * datm,double * dobs,int32 helflag,double mag,double azi_obj,double azi_sun,double azi_moon,double alt_moon,double * dret,char * serr)1695 int32 CALL_CONV swe_heliacal_angle(double tjdut, double *dgeo, double *datm, double *dobs, int32 helflag, double mag, double azi_obj, double azi_sun, double azi_moon, double alt_moon, double *dret, char *serr)
1696 {
1697 if (dgeo[2] < SEI_ECL_GEOALT_MIN || dgeo[2] > SEI_ECL_GEOALT_MAX) {
1698 if (serr != NULL)
1699 sprintf(serr, "location for heliacal events must be between %.0f and %.0f m above sea", SEI_ECL_GEOALT_MIN, SEI_ECL_GEOALT_MAX);
1700 return ERR;
1701 }
1702 swi_set_tid_acc(tjdut, helflag, 0, serr);
1703 default_heliacal_parameters(datm, dgeo, dobs, helflag);
1704 return HeliacalAngle(mag, dobs, azi_obj, alt_moon, azi_moon, tjdut, azi_sun, dgeo, datm, helflag, dret, serr);
1705 }
1706
1707 /*###################################################################
1708 ' AltO [deg]
1709 ' AziO [deg]
1710 ' AltS [deg]
1711 ' AziS [deg]
1712 ' parallax [deg]
1713 ' WidthMoon [deg]
1714 */
WidthMoon(double AltO,double AziO,double AltS,double AziS,double parallax)1715 static double WidthMoon(double AltO, double AziO, double AltS, double AziS, double parallax)
1716 {
1717 /* Yallop 1998, page 3*/
1718 double GeoAltO = AltO + parallax;
1719 return 0.27245 * parallax * (1 + sin(GeoAltO * DEGTORAD) * sin(parallax * DEGTORAD)) * (1 - cos((AltS - GeoAltO) * DEGTORAD) * cos((AziS - AziO) * DEGTORAD));
1720 }
1721
1722 /*###################################################################
1723 ' W [deg]
1724 ' LengthMoon [deg]
1725 */
LengthMoon(double W,double Diamoon)1726 static double LengthMoon(double W, double Diamoon)
1727 {
1728 double Wi, D;
1729 if (Diamoon == 0) Diamoon = AvgRadiusMoon * 2;
1730 Wi = W * 60;
1731 D = Diamoon * 60;
1732 /* Crescent length according: http://calendar.ut.ac.ir/Fa/Crescent/Data/Sultan2005.pdf*/
1733 return (D - 0.3 * (D + Wi) / 2.0 / Wi) / 60.0;
1734 }
1735
1736 /*###################################################################
1737 ' W [deg]
1738 ' GeoARCVact [deg]
1739 ' q [-]
1740 */
qYallop(double W,double GeoARCVact)1741 static double qYallop(double W, double GeoARCVact)
1742 {
1743 double Wi = W * 60;
1744 return (GeoARCVact - (11.8371 - 6.3226 * Wi + 0.7319 * Wi * Wi - 0.1018 * Wi * Wi * Wi)) / 10;
1745 }
1746
1747 /*###################################################################
1748 'A (0,p)
1749 'B (1,q)
1750 'C (0,r)
1751 'D (1,s)
1752 */
crossing(double A,double B,double C,double D)1753 static double crossing(double A, double B, double C, double D)
1754 {
1755 return (C - A) / ((B - A) - (D - C));
1756 }
1757
1758 /*###################################################################*/
DeterTAV(double * dobs,double JDNDaysUT,double * dgeo,double * datm,char * ObjectName,int32 helflag,double * dret,char * serr)1759 static int32 DeterTAV(double *dobs, double JDNDaysUT, double *dgeo, double *datm, char *ObjectName, int32 helflag, double *dret, char *serr)
1760 {
1761 double Magn, AltO, AziS, AziO, AziM, AltM;
1762 double sunra = SunRA(JDNDaysUT, helflag, serr);
1763 if (Magnitude(JDNDaysUT, dgeo, ObjectName, helflag, &Magn, serr) == ERR)
1764 return ERR;
1765 if (ObjectLoc(JDNDaysUT, dgeo, datm, ObjectName, 0, helflag, &AltO, serr) == ERR)
1766 return ERR;
1767 if (ObjectLoc(JDNDaysUT, dgeo, datm, ObjectName, 1, helflag, &AziO, serr) == ERR)
1768 return ERR;
1769 if (strncmp(ObjectName, "moon", 4) == 0) {
1770 AltM = -90;
1771 AziM = 0;
1772 } else {
1773 if (ObjectLoc(JDNDaysUT, dgeo, datm, "moon", 0, helflag, &AltM, serr) == ERR)
1774 return ERR;
1775 if (ObjectLoc(JDNDaysUT, dgeo, datm, "moon", 1, helflag, &AziM, serr) == ERR)
1776 return ERR;
1777 }
1778 if (ObjectLoc(JDNDaysUT, dgeo, datm, "sun", 1, helflag, &AziS, serr) == ERR)
1779 return ERR;
1780 if (TopoArcVisionis(Magn, dobs, AltO, AziO, AltM, AziM, JDNDaysUT, AziS, sunra, dgeo[1], dgeo[2], datm, helflag, dret, serr) == ERR)
1781 return ERR;
1782 return OK;
1783 }
1784
1785 /*###################################################################
1786 ' A y-value at x=1
1787 ' B y-value at x=0
1788 ' C y-value at x=-1
1789 ' x2min minimum for the quadratic function
1790 */
x2min(double A,double B,double C)1791 static double x2min(double A, double B, double C)
1792 {
1793 double term = A + C - 2 * B;
1794 if (term == 0)
1795 return 0;
1796 return -(A - C) / 2.0 / term;
1797 }
1798
1799
1800 /*###################################################################
1801 ' A y-value at x=1
1802 ' B y-value at x=0
1803 ' C y-value at x=-1
1804 ' x
1805 ' y is y-value of quadratic function
1806 */
funct2(double A,double B,double C,double x)1807 static double funct2(double A, double B, double C, double x)
1808 {
1809 return (A + C - 2 * B) / 2.0 * x * x + (A - C) / 2.0 * x + B;
1810 }
1811
strcpy_VBsafe(char * sout,char * sin)1812 static void strcpy_VBsafe(char *sout, char *sin)
1813 {
1814 char *sp, *sp2;
1815 int iw = 0;
1816 sp = sin;
1817 sp2 = sout;
1818 /* note, star name may begin with comma, such as ",zePsc" */
1819 while((isalnum(*sp) || *sp == ' ' || *sp == '-' || *sp == ',') && iw < 30) {
1820 *sp2 = *sp;
1821 sp++; sp2++; iw++;
1822 }
1823 *sp2 = '\0';
1824 }
1825
1826 /*###################################################################
1827 ' JDNDaysUT [JDN]
1828 ' for information about input parameters, see function swe_heliacal_ut().
1829 '
1830 ' output values:
1831 '0=AltO [deg] topocentric altitude of object (unrefracted)
1832 '1=AppAltO [deg] apparent altitude of object (refracted)
1833 '2=GeoAltO [deg] geocentric altitude of object
1834 '3=AziO [deg] azimuth of object
1835 '4=AltS [deg] topocentric altitude of Sun
1836 '5=AziS [deg] azimuth of Sun
1837 '6=TAVact [deg] actual topocentric arcus visionis
1838 '7=ARCVact [deg] actual (geocentric) arcus visionis
1839 '8=DAZact [deg] actual difference between object's and sun's azimuth
1840 '9=ARCLact [deg] actual longitude difference between object and sun
1841 '10=kact [-] extinction coefficient
1842 '11=minTAV [deg] smallest topocentric arcus visionis
1843 '12=TfistVR [JDN] first time object is visible, according to VR
1844 '13=TbVR [JDN] optimum time the object is visible, according to VR
1845 '14=TlastVR [JDN] last time object is visible, according to VR
1846 '15=TbYallop[JDN] best time the object is visible, according to Yallop
1847 '16=WMoon [deg] cresent width of moon
1848 '17=qYal [-] q-test value of Yallop
1849 '18=qCrit [-] q-test criterion of Yallop
1850 '19=ParO [deg] parallax of object
1851 '20 Magn [-] magnitude of object
1852 '21=RiseO [JDN] rise/set time of object
1853 '22=RiseS [JDN] rise/set time of sun
1854 '23=Lag [JDN] rise/set time of object minus rise/set time of sun
1855 '24=TvisVR [JDN] visibility duration
1856 '25=LMoon [deg] cresent length of moon
1857 '26=CVAact [deg]
1858 '27=Illum [%] 'new
1859 '28=CVAact [deg] 'new
1860 '29=MSk [-]
1861 */
swe_heliacal_pheno_ut(double JDNDaysUT,double * dgeo,double * datm,double * dobs,char * ObjectNameIn,int32 TypeEvent,int32 helflag,double * darr,char * serr)1862 int32 CALL_CONV swe_heliacal_pheno_ut(double JDNDaysUT, double *dgeo, double *datm, double *dobs, char *ObjectNameIn, int32 TypeEvent, int32 helflag, double *darr, char *serr)
1863 {
1864 double AziS, AltS, AltS2, AziO, AltO, AltO2, GeoAltO, AppAltO, DAZact, TAVact, ParO, MagnO;
1865 double ARCVact, ARCLact, kact, WMoon, LMoon = 0, qYal, qCrit;
1866 double RiseSetO, RiseSetS, Lag, TbYallop, TfirstVR, TlastVR, TbVR;
1867 double MinTAV = 0, MinTAVact, Ta, Tc, TimeStep, TimePointer, MinTAVoud = 0, DeltaAltoud = 0, DeltaAlt, TvisVR, crosspoint;
1868 double OldestMinTAV, extrax, illum;
1869 double elong, attr[30];
1870 double TimeCheck, LocalminCheck;
1871 int32 retval = OK, RS, Planet;
1872 AS_BOOL noriseO = FALSE;
1873 char ObjectName[AS_MAXCH];
1874 double sunra;
1875 int32 iflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
1876 if (dgeo[2] < SEI_ECL_GEOALT_MIN || dgeo[2] > SEI_ECL_GEOALT_MAX) {
1877 if (serr != NULL)
1878 sprintf(serr, "location for heliacal events must be between %.0f and %.0f m above sea", SEI_ECL_GEOALT_MIN, SEI_ECL_GEOALT_MAX);
1879 return ERR;
1880 }
1881 swi_set_tid_acc(JDNDaysUT, helflag, 0, serr);
1882 sunra = SunRA(JDNDaysUT, helflag, serr);
1883 /* note, the fixed stars functions rewrite the star name. The input string
1884 may be too short, so we have to make sure we have enough space */
1885 strcpy_VBsafe(ObjectName, ObjectNameIn);
1886 tolower_string_star(ObjectName);
1887 default_heliacal_parameters(datm, dgeo, dobs, helflag);
1888 swe_set_topo(dgeo[0], dgeo[1], dgeo[2]);
1889 retval = ObjectLoc(JDNDaysUT, dgeo, datm, "sun", 1, helflag, &AziS, serr);
1890 if (retval == OK)
1891 retval = ObjectLoc(JDNDaysUT, dgeo, datm, "sun", 0, helflag, &AltS, serr);
1892 if (retval == OK)
1893 retval = ObjectLoc(JDNDaysUT, dgeo, datm, ObjectName, 1, helflag, &AziO, serr);
1894 if (retval == OK)
1895 retval = ObjectLoc(JDNDaysUT, dgeo, datm, ObjectName, 0, helflag, &AltO, serr);
1896 if (retval == OK)
1897 retval = ObjectLoc(JDNDaysUT, dgeo, datm, ObjectName, 7, helflag, &GeoAltO, serr);
1898 if (retval == ERR)
1899 return ERR;
1900 AppAltO = AppAltfromTopoAlt(AltO, datm[1], datm[0], helflag);
1901 DAZact = AziS - AziO;
1902 TAVact = AltO - AltS;
1903 /*this parallax seems to be somewhat smaller then in Yallop and SkyMap! Needs to be studied*/
1904 ParO = GeoAltO - AltO;
1905 if (Magnitude(JDNDaysUT, dgeo, ObjectName, helflag, &MagnO, serr) == ERR)
1906 return ERR;
1907 ARCVact = TAVact + ParO;
1908 ARCLact = acos(cos(ARCVact * DEGTORAD) * cos(DAZact * DEGTORAD)) / DEGTORAD;
1909 Planet = DeterObject(ObjectName);
1910 if (Planet == -1) {
1911 elong = ARCLact;
1912 illum = 100;
1913 } else {
1914 retval = swe_pheno_ut(JDNDaysUT, Planet, iflag|(SEFLG_TOPOCTR|SEFLG_EQUATORIAL), attr, serr);
1915 if (retval == ERR) return ERR;
1916 elong = attr[2];
1917 illum = attr[1] * 100;
1918 }
1919 kact = kt(AltS, sunra, dgeo[1], dgeo[2], datm[1], datm[2], datm[3], 4, serr);
1920 if ((0)) {
1921 darr[26] = kR(AltS, dgeo[2]);
1922 darr[27] = kW(dgeo[2], datm[1], datm[2]);
1923 darr[28] = kOZ(AltS, sunra, dgeo[1]);
1924 darr[29] = ka(AltS, sunra, dgeo[1], dgeo[2], datm[1], datm[2], datm[3], serr);
1925 darr[30] = darr[26] + darr[27] + darr[28] + darr[29];
1926 }
1927 WMoon = 0;
1928 qYal = 0;
1929 qCrit = 0;
1930 LMoon = 0;
1931 if (Planet == SE_MOON) {
1932 WMoon = WidthMoon(AltO, AziO, AltS, AziS, ParO);
1933 LMoon = LengthMoon(WMoon, 0);
1934 qYal = qYallop(WMoon, ARCVact);
1935 if (qYal > 0.216) qCrit = 1; /* A */
1936 if (qYal < 0.216 && qYal > -0.014) qCrit = 2; /* B */
1937 if (qYal < -0.014 && qYal > -0.16) qCrit = 3; /* C */
1938 if (qYal < -0.16 && qYal > -0.232) qCrit = 4; /* D */
1939 if (qYal < -0.232 && qYal > -0.293) qCrit = 5; /* E */
1940 if (qYal < -0.293) qCrit = 6; /* F */
1941 }
1942 /*determine if rise or set event*/
1943 RS = 2;
1944 if (TypeEvent == 1 || TypeEvent == 4) RS = 1;
1945 retval = RiseSet(JDNDaysUT - 4.0 / 24.0, dgeo, datm, "sun", RS, helflag, 0, &RiseSetS, serr);
1946 if (retval == ERR)
1947 return ERR;
1948 retval = RiseSet(JDNDaysUT - 4.0 / 24.0, dgeo, datm, ObjectName, RS, helflag, 0, &RiseSetO, serr);
1949 if (retval == ERR)
1950 return ERR;
1951 TbYallop = TJD_INVALID;
1952 if (retval == -2) { /* object does not rise or set */
1953 Lag = 0;
1954 noriseO = TRUE;
1955 } else {
1956 Lag = RiseSetO - RiseSetS;
1957 if (Planet == SE_MOON)
1958 TbYallop = (RiseSetO * 4 + RiseSetS * 5) / 9.0;
1959 }
1960 if ((TypeEvent == 3 || TypeEvent == 4) && (Planet == -1 || Planet >= SE_MARS)) {
1961 TfirstVR = TJD_INVALID;
1962 TbVR = TJD_INVALID;
1963 TlastVR = TJD_INVALID;
1964 TvisVR = 0;
1965 MinTAV = 0;
1966 goto output_heliacal_pheno;
1967 }
1968 /* If HPheno >= 11 And HPheno <= 14 Or HPheno = 24 Then*/
1969 /*te bepalen m.b.v. walkthrough*/
1970 MinTAVact = 199;
1971 DeltaAlt = 0;
1972 OldestMinTAV = 0;
1973 Ta = 0;
1974 Tc = 0;
1975 TbVR = 0;
1976 TimeStep = -TimeStepDefault / 24.0 / 60.0;
1977 if (RS == 2) TimeStep = -TimeStep;
1978 TimePointer = RiseSetS - TimeStep;
1979 do {
1980 TimePointer = TimePointer + TimeStep;
1981 OldestMinTAV = MinTAVoud;
1982 MinTAVoud = MinTAVact;
1983 DeltaAltoud = DeltaAlt;
1984 retval = ObjectLoc(TimePointer, dgeo, datm, "sun", 0, helflag, &AltS2, serr);
1985 if (retval == OK)
1986 retval = ObjectLoc(TimePointer, dgeo, datm, ObjectName, 0, helflag, &AltO2, serr);
1987 if (retval != OK)
1988 return ERR;
1989 DeltaAlt = AltO2 - AltS2;
1990 if (DeterTAV(dobs, TimePointer, dgeo, datm, ObjectName, helflag, &MinTAVact, serr) == ERR)
1991 return ERR;
1992 if (MinTAVoud < MinTAVact && TbVR == 0) {
1993 /* determine if this is a local minimum with object still above horizon*/
1994 TimeCheck = TimePointer + Sgn(TimeStep) * LocalMinStep / 24.0 / 60.0;
1995 if (RiseSetO != 0) {
1996 if (TimeStep > 0)
1997 TimeCheck = mymin(TimeCheck, RiseSetO);
1998 else
1999 TimeCheck = mymax(TimeCheck, RiseSetO);
2000 }
2001 if (DeterTAV(dobs, TimeCheck, dgeo, datm, ObjectName, helflag, &LocalminCheck, serr) == ERR)
2002 return ERR;
2003 if (LocalminCheck > MinTAVact) {
2004 extrax = x2min(MinTAVact, MinTAVoud, OldestMinTAV);
2005 TbVR = TimePointer - (1 - extrax) * TimeStep;
2006 MinTAV = funct2(MinTAVact, MinTAVoud, OldestMinTAV, extrax);
2007 }
2008 }
2009 if (DeltaAlt > MinTAVact && Tc == 0 && TbVR == 0) {
2010 crosspoint = crossing(DeltaAltoud, DeltaAlt, MinTAVoud, MinTAVact);
2011 Tc = TimePointer - TimeStep * (1 - crosspoint);
2012 }
2013 if (DeltaAlt < MinTAVact && Ta == 0 && Tc != 0) {
2014 crosspoint = crossing(DeltaAltoud, DeltaAlt, MinTAVoud, MinTAVact);
2015 Ta = TimePointer - TimeStep * (1 - crosspoint);
2016 }
2017 } while (fabs(TimePointer - RiseSetS) <= MaxTryHours / 24.0 && Ta == 0 && !((TbVR != 0 && (TypeEvent == 3 || TypeEvent == 4) && (strncmp(ObjectName, "moon", 4) != 0 && strncmp(ObjectName, "venus", 5) != 0 && strncmp(ObjectName, "mercury", 7) != 0))));
2018 if (RS == 2) {
2019 TfirstVR = Tc;
2020 TlastVR = Ta;
2021 } else {
2022 TfirstVR = Ta;
2023 TlastVR = Tc;
2024 }
2025 if (TfirstVR == 0 && TlastVR == 0) {
2026 if (RS == 1)
2027 TfirstVR = TbVR - 0.000001;
2028 else
2029 TlastVR = TbVR + 0.000001;
2030 }
2031 if (!noriseO) {
2032 if (RS == 1)
2033 TfirstVR = mymax(TfirstVR, RiseSetO);
2034 else
2035 TlastVR = mymin(TlastVR, RiseSetO);
2036 }
2037 TvisVR = TJD_INVALID; /*"#NA!" */
2038 if (TlastVR != 0 && TfirstVR != 0)
2039 TvisVR = TlastVR - TfirstVR;
2040 if (TlastVR == 0) TlastVR = TJD_INVALID; /*"#NA!" */
2041 if (TbVR == 0) TbVR = TJD_INVALID; /*"#NA!" */
2042 if (TfirstVR == 0) TfirstVR = TJD_INVALID; /*"#NA!" */
2043 output_heliacal_pheno:
2044 /*End If*/
2045 darr[0] = AltO;
2046 darr[1] = AppAltO;
2047 darr[2] = GeoAltO;
2048 darr[3] = AziO;
2049 darr[4] = AltS;
2050 darr[5] = AziS;
2051 darr[6] = TAVact;
2052 darr[7] = ARCVact;
2053 darr[8] = DAZact;
2054 darr[9] = ARCLact;
2055 darr[10] = kact;
2056 darr[11] = MinTAV;
2057 darr[12] = TfirstVR;
2058 darr[13] = TbVR;
2059 darr[14] = TlastVR;
2060 darr[15] = TbYallop;
2061 darr[16] = WMoon;
2062 darr[17] = qYal;
2063 darr[18] = qCrit;
2064 darr[19] = ParO;
2065 darr[20] = MagnO;
2066 darr[21] = RiseSetO;
2067 darr[22] = RiseSetS;
2068 darr[23] = Lag;
2069 darr[24] = TvisVR;
2070 darr[25] = LMoon;
2071 darr[26] = elong;
2072 darr[27] = illum;
2073 return OK;
2074 }
2075
2076 #if 0
2077 int32 HeliacalJDut(double JDNDaysUTStart, double Age, double SN, double Lat, double Longitude, double HeightEye, double Temperature, double Pressure, double RH, double VR, char *ObjectName, int TypeEvent, char *AVkind, double *dret, char *serr)
2078 {
2079 double dgeo[3], datm[4], dobs[6];
2080 int32 helflag = SE_HELFLAG_HIGH_PRECISION;
2081 helflag |= SE_HELFLAG_AVKIND_VR;
2082 dgeo[0] = Longitude;
2083 dgeo[1] = Lat;
2084 dgeo[2] = HeightEye;
2085 datm[0] = Pressure;
2086 datm[1] = Temperature;
2087 datm[2] = RH;
2088 datm[3] = VR;
2089 dobs[0] = Age;
2090 dobs[1] = SN;
2091 return swe_heliacal_ut(JDNDaysUTStart, dgeo, datm, dobs, ObjectName, TypeEvent, helflag, dret, serr);
2092 }
2093 #endif
2094
get_synodic_period(int Planet)2095 static double get_synodic_period(int Planet)
2096 {
2097 /* synodic periods from:
2098 * Kelley/Milone/Aveni, "Exploring ancient Skies", p. 43. */
2099 switch(Planet) {
2100 case SE_MOON: return 29.530588853;
2101 case SE_MERCURY: return 115.8775;
2102 case SE_VENUS: return 583.9214;
2103 case SE_MARS: return 779.9361;
2104 case SE_JUPITER: return 398.8840;
2105 case SE_SATURN: return 378.0919;
2106 case SE_URANUS: return 369.6560;
2107 case SE_NEPTUNE: return 367.4867;
2108 case SE_PLUTO: return 366.7207;
2109 }
2110 return 366; /* for stars and default for far away planets */
2111 }
2112
2113 /*###################################################################*/
moon_event_arc_vis(double JDNDaysUTStart,double * dgeo,double * datm,double * dobs,int32 TypeEvent,int32 helflag,double * dret,char * serr)2114 static int32 moon_event_arc_vis(double JDNDaysUTStart, double *dgeo, double *datm, double *dobs, int32 TypeEvent, int32 helflag, double *dret, char *serr)
2115 {
2116 double x[20], MinTAV, MinTAVoud, OldestMinTAV;
2117 double phase1, phase2, JDNDaysUT, JDNDaysUTi;
2118 double tjd_moonevent, tjd_moonevent_start;
2119 double DeltaAltoud, TimeCheck, LocalminCheck;
2120 double AltS, AltO, DeltaAlt = 90;
2121 char ObjectName[30];
2122 int32 iflag, Daystep, goingup, Planet, retval;
2123 int32 avkind = helflag & SE_HELFLAG_AVKIND;
2124 int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
2125 dret[0] = JDNDaysUTStart; /* will be returned in error case */
2126 if (avkind == 0)
2127 avkind = SE_HELFLAG_AVKIND_VR;
2128 if (avkind != SE_HELFLAG_AVKIND_VR) {
2129 if (serr != NULL)
2130 strcpy(serr, "error: in valid AV kind for the moon");
2131 return ERR;
2132 }
2133 if (TypeEvent == 1 || TypeEvent == 2) {
2134 if (serr != NULL)
2135 strcpy(serr, "error: the moon has no morning first or evening last");
2136 return ERR;
2137 }
2138 strcpy(ObjectName, "moon");
2139 Planet = SE_MOON;
2140 iflag = SEFLG_TOPOCTR | SEFLG_EQUATORIAL | epheflag;
2141 if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
2142 iflag |= SEFLG_NONUT|SEFLG_TRUEPOS;
2143 Daystep = 1;
2144 if (TypeEvent == 3) {
2145 /*morning last */
2146 TypeEvent = 2;
2147 } else {
2148 /*evening first*/
2149 TypeEvent = 1;
2150 Daystep = -Daystep;
2151 }
2152 /* check Synodic/phase Period */
2153 JDNDaysUT = JDNDaysUTStart;
2154 /* start 30 days later if TypeEvent=4 (1) */
2155 if (TypeEvent == 1) JDNDaysUT = JDNDaysUT + 30;
2156 /* determination of new moon date */
2157 swe_pheno_ut(JDNDaysUT, Planet, iflag, x, serr);
2158 phase2 = x[0];
2159 goingup = 0;
2160 do {
2161 JDNDaysUT = JDNDaysUT + Daystep;
2162 phase1 = phase2;
2163 swe_pheno_ut(JDNDaysUT, Planet, iflag, x, serr);
2164 phase2 = x[0];
2165 if (phase2 > phase1)
2166 goingup = 1;
2167 } while (goingup == 0 || (goingup == 1 && (phase2 > phase1)));
2168 /* fix the date to get the day with the smallest phase (nwest moon) */
2169 JDNDaysUT = JDNDaysUT - Daystep;
2170 /* initialize the date to look for set */
2171 JDNDaysUTi = JDNDaysUT;
2172 JDNDaysUT = JDNDaysUT - Daystep;
2173 MinTAVoud = 199;
2174 do {
2175 JDNDaysUT = JDNDaysUT + Daystep;
2176 if ((retval = RiseSet(JDNDaysUT, dgeo, datm, ObjectName, TypeEvent, helflag, 0, &tjd_moonevent, serr)) != OK)
2177 return retval;
2178 tjd_moonevent_start = tjd_moonevent;
2179 MinTAV = 199;
2180 OldestMinTAV = MinTAV;
2181 do {
2182 OldestMinTAV = MinTAVoud;
2183 MinTAVoud = MinTAV;
2184 DeltaAltoud = DeltaAlt;
2185 tjd_moonevent = tjd_moonevent - 1.0 / 60.0 / 24.0 * Sgn(Daystep);
2186 if (ObjectLoc(tjd_moonevent, dgeo, datm, "sun", 0, helflag, &AltS, serr) == ERR)
2187 return ERR;
2188 if (ObjectLoc(tjd_moonevent, dgeo, datm, ObjectName, 0, helflag, &AltO, serr) == ERR)
2189 return ERR;
2190 DeltaAlt = AltO - AltS;
2191 if (DeterTAV(dobs, tjd_moonevent, dgeo, datm, ObjectName, helflag, &MinTAV, serr) == ERR)
2192 return ERR;
2193 TimeCheck = tjd_moonevent - LocalMinStep / 60.0 / 24.0 * Sgn(Daystep);
2194 if (DeterTAV(dobs, TimeCheck, dgeo, datm, ObjectName, helflag, &LocalminCheck, serr) == ERR)
2195 return ERR;
2196 /*printf("%f, %f <= %f\n", tjd_moonevent, MinTAV, MinTAVoud);*/
2197 /* while (MinTAV <= MinTAVoud && fabs(tjd_moonevent - tjd_moonevent_start) < 120.0 / 60.0 / 24.0);*/
2198 } while ((MinTAV <= MinTAVoud || LocalminCheck < MinTAV) && fabs(tjd_moonevent - tjd_moonevent_start) < 120.0 / 60.0 / 24.0);
2199 /* while (DeltaAlt < MinTAVoud && fabs(JDNDaysUT - JDNDaysUTi) < 15);*/
2200 } while (DeltaAltoud < MinTAVoud && fabs(JDNDaysUT - JDNDaysUTi) < 15);
2201 if (fabs(JDNDaysUT - JDNDaysUTi) < 15) {
2202 tjd_moonevent += (1 - x2min(MinTAV, MinTAVoud, OldestMinTAV)) * Sgn(Daystep) / 60.0 / 24.0;
2203 } else {
2204 strcpy(serr, "no date found for lunar event");
2205 return ERR;
2206 }
2207 dret[0] = tjd_moonevent;
2208 return OK;
2209 }
2210
heliacal_ut_arc_vis(double JDNDaysUTStart,double * dgeo,double * datm,double * dobs,char * ObjectName,int32 TypeEventIn,int32 helflag,double * dret,char * serr_ret)2211 static int32 heliacal_ut_arc_vis(double JDNDaysUTStart, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 TypeEventIn, int32 helflag, double *dret, char *serr_ret)
2212 {
2213 double x[6];
2214 double xin[2];
2215 double xaz[2];
2216 double dang[3];
2217 double objectmagn = 0, maxlength, DayStep;
2218 double JDNDaysUT, JDNDaysUTfinal, JDNDaysUTstep, JDNDaysUTstepoud, JDNarcvisUT, tjd_tt, tret, OudeDatum, JDNDaysUTinp = JDNDaysUTStart, JDNDaysUTtijd;
2219 double ArcusVis, ArcusVisDelta, ArcusVisPto, ArcusVisDeltaoud;
2220 double Trise, sunsangle, Theliacal, Tdelta, Angle;
2221 double TimeStep, TimePointer, OldestMinTAV, MinTAVoud, MinTAVact, extrax, TbVR = 0;
2222 double AziS, AltS, AziO, AltO, DeltaAlt;
2223 double direct, Pressure, Temperature, d;
2224 int32 epheflag, retval = OK;
2225 int32 iflag, Planet, eventtype;
2226 int32 TypeEvent = TypeEventIn;
2227 int doneoneday;
2228 char serr[AS_MAXCH];
2229 *dret = JDNDaysUTStart;
2230 *serr = '\0';
2231 Planet = DeterObject(ObjectName);
2232 Pressure = datm[0];
2233 Temperature = datm[1];
2234 /* determine Magnitude of star*/
2235 if ((retval = Magnitude(JDNDaysUTStart, dgeo, ObjectName, helflag, &objectmagn, serr)) == ERR)
2236 goto swe_heliacal_err;
2237 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
2238 iflag = SEFLG_TOPOCTR | SEFLG_EQUATORIAL | epheflag;
2239 if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
2240 iflag |= SEFLG_NONUT | SEFLG_TRUEPOS;
2241 /* start values for search of heliacal rise
2242 * maxlength = phase period in days, smaller than minimal synodic period */
2243 /* days per step (for heliacal rise) in power of two */
2244 switch(Planet) {
2245 case SE_MERCURY:
2246 DayStep = 1; maxlength = 100; break;
2247 case SE_VENUS:
2248 DayStep = 64; maxlength = 384; break;
2249 case SE_MARS:
2250 DayStep = 128; maxlength = 640; break;
2251 case SE_JUPITER:
2252 DayStep = 64; maxlength = 384; break;
2253 case SE_SATURN:
2254 DayStep = 64; maxlength = 256; break;
2255 default:
2256 DayStep = 64; maxlength = 256; break;
2257 }
2258 /* heliacal setting */
2259 eventtype = TypeEvent;
2260 if (eventtype == 2) DayStep = -DayStep;
2261 /* acronychal setting */
2262 if (eventtype == 4) {
2263 eventtype = 1;
2264 DayStep = -DayStep;
2265 }
2266 /* acronychal rising */
2267 if (eventtype == 3) eventtype = 2;
2268 eventtype |= SE_BIT_DISC_CENTER;
2269 /* normalize the maxlength to the step size */
2270 {
2271 /* check each Synodic/phase Period */
2272 JDNDaysUT = JDNDaysUTStart;
2273 /* make sure one can find an event on the just after the JDNDaysUTStart */
2274 JDNDaysUTfinal = JDNDaysUT + maxlength;
2275 JDNDaysUT = JDNDaysUT - 1;
2276 if (DayStep < 0) {
2277 JDNDaysUTtijd = JDNDaysUT;
2278 JDNDaysUT = JDNDaysUTfinal;
2279 JDNDaysUTfinal = JDNDaysUTtijd;
2280 }
2281 /* prepair the search */
2282 JDNDaysUTstep = JDNDaysUT - DayStep;
2283 doneoneday = 0;
2284 ArcusVisDelta = 199;
2285 ArcusVisPto = -5.55;
2286 do { /* this is a do {} while() loop */
2287 if (fabs(DayStep) == 1) doneoneday = 1;
2288 do { /* this is a do {} while() loop */
2289 /* init search for heliacal rise */
2290 JDNDaysUTstepoud = JDNDaysUTstep;
2291 ArcusVisDeltaoud = ArcusVisDelta;
2292 JDNDaysUTstep = JDNDaysUTstep + DayStep;
2293 /* determine rise/set time */
2294 if ((retval = my_rise_trans(JDNDaysUTstep, SE_SUN, "", eventtype, helflag, dgeo, datm, &tret, serr)) == ERR)
2295 goto swe_heliacal_err;
2296 /* determine time compensation to get Sun's altitude at heliacal rise */
2297 tjd_tt = tret + swe_deltat_ex(tret, epheflag, serr);
2298 if ((retval = swe_calc(tjd_tt, SE_SUN, iflag, x, serr)) == ERR)
2299 goto swe_heliacal_err;
2300 xin[0] = x[0];
2301 xin[1] = x[1];
2302 swe_azalt(tret, SE_EQU2HOR, dgeo, Pressure, Temperature, xin, xaz);
2303 Trise = HourAngle(xaz[1], x[1], dgeo[1]);
2304 sunsangle = ArcusVisPto;
2305 if (helflag & SE_HELFLAG_AVKIND_MIN7) sunsangle = -7;
2306 if (helflag & SE_HELFLAG_AVKIND_MIN9) sunsangle = -9;
2307 Theliacal = HourAngle(sunsangle, x[1], dgeo[1]);
2308 Tdelta = Theliacal - Trise;
2309 if (TypeEvent == 2 || TypeEvent== 3) Tdelta = -Tdelta;
2310 /* determine appr.time when sun is at the wanted Sun's altitude */
2311 JDNarcvisUT = tret - Tdelta / 24;
2312 tjd_tt = JDNarcvisUT + swe_deltat_ex(JDNarcvisUT, epheflag, serr);
2313 /* determine Sun's position */
2314 if ((retval = swe_calc(tjd_tt, SE_SUN, iflag, x, serr)) == ERR)
2315 goto swe_heliacal_err;
2316 xin[0] = x[0];
2317 xin[1] = x[1];
2318 swe_azalt(JDNarcvisUT, SE_EQU2HOR, dgeo, Pressure, Temperature, xin, xaz);
2319 AziS = xaz[0] + 180;
2320 if (AziS >= 360) AziS = AziS - 360;
2321 AltS = xaz[1];
2322 /* determine Moon's position */
2323 #if 0
2324 double AltM, AziM;
2325 if ((retval = swe_calc(tjd_tt, SE_MOON, iflag, x, serr)) == ERR)
2326 goto swe_heliacal_err;
2327 xin[0] = x[0];
2328 xin[1] = x[1];
2329 swe_azalt(JDNarcvisUT, SE_EQU2HOR, dgeo, Pressure, Temperature, xin, xaz);
2330 AziM = xaz[0] + 180;
2331 if (AziM >= 360) AziM = AziM - 360;
2332 AltM = xaz[1];
2333 #endif
2334 /* determine object's position */
2335 if (Planet != -1) {
2336 if ((retval = swe_calc(tjd_tt, Planet, iflag, x, serr)) == ERR)
2337 goto swe_heliacal_err;
2338 /* determine magnitude of Planet */
2339 if ((retval = Magnitude(JDNarcvisUT, dgeo, ObjectName, helflag, &objectmagn, serr)) == ERR)
2340 goto swe_heliacal_err;
2341 } else {
2342 if ((retval = call_swe_fixstar(ObjectName, tjd_tt, iflag, x, serr)) == ERR)
2343 goto swe_heliacal_err;
2344 }
2345 xin[0] = x[0];
2346 xin[1] = x[1];
2347 swe_azalt(JDNarcvisUT, SE_EQU2HOR, dgeo, Pressure, Temperature, xin, xaz);
2348 AziO = xaz[0] + 180;
2349 if (AziO >= 360) AziO = AziO - 360;
2350 AltO = xaz[1];
2351 /* determine arcusvisionis */
2352 DeltaAlt = AltO - AltS;
2353 /*if ((retval = HeliacalAngle(objectmagn, dobs, AziO, AltM, AziM, JDNarcvisUT, AziS, dgeo, datm, helflag, dang, serr)) == ERR)*/
2354 if ((retval = HeliacalAngle(objectmagn, dobs, AziO, -1, 0, JDNarcvisUT, AziS, dgeo, datm, helflag, dang, serr)) == ERR)
2355 goto swe_heliacal_err;
2356 ArcusVis = dang[1];
2357 ArcusVisPto = dang[2];
2358 ArcusVisDelta = DeltaAlt - ArcusVis;
2359 /*} while (((ArcusVisDeltaoud > 0 && ArcusVisDelta < 0) || ArcusVisDelta < 0) && (JDNDaysUTfinal - JDNDaysUTstep) * Sgn(DayStep) > 0);*/
2360 } while ((ArcusVisDeltaoud > 0 || ArcusVisDelta < 0) && (JDNDaysUTfinal - JDNDaysUTstep) * Sgn(DayStep) > 0);
2361 if (doneoneday == 0 && (JDNDaysUTfinal - JDNDaysUTstep) * Sgn(DayStep) > 0) {
2362 /* go back to date before heliacal altitude */
2363 ArcusVisDelta = ArcusVisDeltaoud;
2364 DayStep = ((int) (fabs(DayStep) / 2.0)) * Sgn(DayStep);
2365 JDNDaysUTstep = JDNDaysUTstepoud;
2366 }
2367 } while (doneoneday == 0 && (JDNDaysUTfinal - JDNDaysUTstep) * Sgn(DayStep) > 0);
2368 }
2369 d = (JDNDaysUTfinal - JDNDaysUTstep) * Sgn(DayStep);
2370 if (d <= 0 || d >= maxlength) {
2371 dret[0] = JDNDaysUTinp; /* no date found, just return input */
2372 retval = -2; /* marks "not found" within synodic period */
2373 sprintf(serr, "heliacal event not found within maxlength %f\n", maxlength);
2374 goto swe_heliacal_err;
2375 }
2376 #if 0
2377 if (helflag & SE_HELFLAG_AVKIND_VR) {
2378 double darr[40];
2379 if (swe_heliacal_pheno_ut(JDNarcvisUT, dgeo, datm, dobs, ObjectName, TypeEvent, helflag, darr, serr) != OK)
2380 return ERR;
2381 JDNarcvisUT = darr[13];
2382 }
2383 }
2384 #endif
2385 direct = TimeStepDefault / 24.0 / 60.0;
2386 if (DayStep < 0) direct = -direct;
2387 if (helflag & SE_HELFLAG_AVKIND_VR) {
2388 /*te bepalen m.b.v. walkthrough*/
2389 TimeStep = direct;
2390 TbVR = 0;
2391 TimePointer = JDNarcvisUT;
2392 if (DeterTAV(dobs, TimePointer, dgeo, datm, ObjectName, helflag, &OldestMinTAV, serr) == ERR)
2393 return ERR;
2394 TimePointer = TimePointer + TimeStep;
2395 if (DeterTAV(dobs, TimePointer, dgeo, datm, ObjectName, helflag, &MinTAVoud, serr) == ERR)
2396 return ERR;
2397 if (MinTAVoud > OldestMinTAV) {
2398 TimePointer = JDNarcvisUT;
2399 TimeStep = -TimeStep;
2400 MinTAVact = OldestMinTAV;
2401 } else {
2402 MinTAVact = MinTAVoud;
2403 MinTAVoud = OldestMinTAV;
2404 }
2405 /*TimePointer = TimePointer - Timestep*/
2406 do {
2407 TimePointer = TimePointer + TimeStep;
2408 OldestMinTAV = MinTAVoud;
2409 MinTAVoud = MinTAVact;
2410 if (DeterTAV(dobs, TimePointer, dgeo, datm, ObjectName, helflag, &MinTAVact, serr) == ERR)
2411 return ERR;
2412 if (MinTAVoud < MinTAVact) {
2413 extrax = x2min(MinTAVact, MinTAVoud, OldestMinTAV);
2414 TbVR = TimePointer - (1 - extrax) * TimeStep;
2415 }
2416 } while (TbVR == 0);
2417 JDNarcvisUT = TbVR;
2418 }
2419 /*if (strncmp(AVkind, "pto", 3) == 0) */
2420 if (helflag & SE_HELFLAG_AVKIND_PTO) {
2421 do {
2422 OudeDatum = JDNarcvisUT;
2423 JDNarcvisUT = JDNarcvisUT - direct;
2424 tjd_tt = JDNarcvisUT + swe_deltat_ex(JDNarcvisUT, epheflag, serr);
2425 if (Planet != -1) {
2426 if ((retval = swe_calc(tjd_tt, Planet, iflag, x, serr)) == ERR)
2427 goto swe_heliacal_err;
2428 } else {
2429 if ((retval = call_swe_fixstar(ObjectName, tjd_tt, iflag, x, serr)) == ERR)
2430 goto swe_heliacal_err;
2431 }
2432 xin[0] = x[0];
2433 xin[1] = x[1];
2434 swe_azalt(JDNarcvisUT, SE_EQU2HOR, dgeo, Pressure, Temperature, xin, xaz);
2435 Angle = xaz[1];
2436 } while (Angle > 0);
2437 JDNarcvisUT = (JDNarcvisUT + OudeDatum) / 2.0;
2438 }
2439 if (JDNarcvisUT < -9999999 || JDNarcvisUT > 9999999) {
2440 dret[0] = JDNDaysUT; /* no date found, just return input */
2441 strcpy(serr, "no heliacal date found");
2442 retval = ERR;
2443 goto swe_heliacal_err;
2444 }
2445 dret[0] = JDNarcvisUT;
2446 swe_heliacal_err:
2447 if (serr_ret != NULL && *serr != '\0')
2448 strcpy(serr_ret, serr);
2449 return retval;
2450 }
2451
get_asc_obl(double tjd,int32 ipl,char * star,int32 iflag,double * dgeo,AS_BOOL desc_obl,double * daop,char * serr)2452 static int32 get_asc_obl(double tjd, int32 ipl, char *star, int32 iflag, double
2453 *dgeo, AS_BOOL desc_obl, double *daop, char *serr)
2454 {
2455 int32 retval;
2456 int32 epheflag = iflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
2457 double x[6], adp;
2458 char s[AS_MAXCH];
2459 char star2[AS_MAXCH];
2460 strcpy(star2, star);
2461 if (ipl == -1) {
2462 if ((retval = swe_fixstar(star2, tjd, epheflag | SEFLG_EQUATORIAL, x, serr)) == ERR)
2463 return ERR;
2464 } else {
2465 if ((retval = swe_calc(tjd, ipl, epheflag | SEFLG_EQUATORIAL, x, serr)) == ERR)
2466 return ERR;
2467 }
2468 adp = tan(dgeo[1] * DEGTORAD) * tan(x[1] * DEGTORAD);
2469 if (fabs(adp) > 1) {
2470 if (star != NULL && *star != '\0')
2471 strcpy(s, star);
2472 else
2473 swe_get_planet_name(ipl, s);
2474 sprintf(serr, "%s is circumpolar, cannot calculate heliacal event", s);
2475 return -2;
2476 }
2477 adp = asin(adp) / DEGTORAD;
2478 if (desc_obl)
2479 *daop = x[0] + adp;
2480 else
2481 *daop = x[0] - adp;
2482 *daop = swe_degnorm(*daop);
2483 return OK;
2484 }
2485
2486 #if 0
2487 static int32 get_asc_obl_old(double tjd, int32 ipl, char *star, int32 iflag, double *dgeo, AS_BOOL desc_obl, double *daop, char *serr)
2488 {
2489 int32 retval;
2490 int32 epheflag = iflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
2491 double x[6], adp;
2492 char s[AS_MAXCH];
2493 if (star != NULL && *star != '\0') {
2494 if ((retval = call_swe_fixstar(star, tjd, epheflag | SEFLG_EQUATORIAL, x, serr)) == ERR)
2495 return ERR;
2496 } else {
2497 if ((retval = swe_calc(tjd, ipl, epheflag | SEFLG_EQUATORIAL, x, serr)) == ERR)
2498 return ERR;
2499 }
2500 adp = tan(dgeo[1] * DEGTORAD) * tan(x[1] * DEGTORAD);
2501 if (fabs(adp) > 1) {
2502 if (star != NULL && *star != '\0')
2503 strcpy(s, star);
2504 else
2505 swe_get_planet_name(ipl, s);
2506 sprintf(serr, "%s is circumpolar, cannot calculate heliacal event", s);
2507 return -2;
2508 }
2509 adp = asin(adp) / DEGTORAD;
2510 if (desc_obl)
2511 *daop = x[0] + adp;
2512 else
2513 *daop = x[0] - adp;
2514 *daop = swe_degnorm(*daop);
2515 return OK;
2516 }
2517 #endif
2518
get_asc_obl_diff(double tjd,int32 ipl,char * star,int32 iflag,double * dgeo,AS_BOOL desc_obl,AS_BOOL is_acronychal,double * dsunpl,char * serr)2519 static int32 get_asc_obl_diff(double tjd, int32 ipl, char *star, int32 iflag, double *dgeo, AS_BOOL desc_obl, AS_BOOL is_acronychal, double *dsunpl, char *serr)
2520 {
2521 int32 retval = OK;
2522 double aosun, aopl;
2523 /* ascensio obliqua of sun */
2524 retval = get_asc_obl(tjd, SE_SUN, "", iflag, dgeo, desc_obl, &aosun, serr);
2525 if (retval != OK)
2526 return retval;
2527 if (is_acronychal) {
2528 if (desc_obl == TRUE)
2529 desc_obl = FALSE;
2530 else
2531 desc_obl = TRUE;
2532 }
2533 /* ascensio obliqua of body */
2534 retval = get_asc_obl(tjd, ipl, star, iflag, dgeo, desc_obl, &aopl, serr);
2535 if (retval != OK)
2536 return retval;
2537 *dsunpl = swe_degnorm(aosun - aopl);
2538 if (is_acronychal)
2539 *dsunpl = swe_degnorm(*dsunpl - 180);
2540 if (*dsunpl > 180) *dsunpl -= 360;
2541 return OK;
2542 }
2543
2544 #if 0
2545 static int32 get_asc_obl_diff_old(double tjd, int32 ipl, char *star, int32 iflag, double *dgeo, AS_BOOL desc_obl, double *dsunpl, char *serr)
2546 {
2547 int32 retval = OK;
2548 double aosun, aopl;
2549 /* ascensio obliqua of sun */
2550 retval = get_asc_obl(tjd, SE_SUN, "", iflag, dgeo, desc_obl, &aosun, serr);
2551 if (retval != OK)
2552 return retval;
2553 /* ascensio obliqua of body */
2554 retval = get_asc_obl(tjd, ipl, star, iflag, dgeo, desc_obl, &aopl, serr);
2555 if (retval != OK)
2556 return retval;
2557 *dsunpl = swe_degnorm(aosun - aopl);
2558 return OK;
2559 }
2560 #endif
2561
2562 /* times of
2563 * - superior and inferior conjunction (Mercury and Venus)
2564 * - conjunction and opposition (ipl >= Mars)
2565 */
2566 static const double tcon[] =
2567 {
2568 0, 0,
2569 2451550, 2451550, /* Moon */
2570 2451604, 2451670, /* Mercury */
2571 2451980, 2452280, /* Venus */
2572 2451727, 2452074, /* Mars */
2573 2451673, 2451877, /* Jupiter */
2574 2451675, 2451868, /* Saturn */
2575 2451581, 2451768, /* Uranus */
2576 2451568, 2451753, /* Neptune */
2577 };
2578
find_conjunct_sun(double tjd_start,int32 ipl,int32 helflag,int32 TypeEvent,double * tjd,char * serr)2579 static int32 find_conjunct_sun(double tjd_start, int32 ipl, int32 helflag, int32 TypeEvent, double *tjd, char *serr)
2580 {
2581 int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
2582 int i;
2583 double tjdcon, tjd0, ds, dsynperiod, x[6], xs[6], daspect = 0;
2584 if (ipl >= SE_MARS && TypeEvent >= 3)
2585 daspect = 180;
2586 i = (TypeEvent - 1) / 2 + ipl * 2;
2587 tjd0 = tcon[i];
2588 dsynperiod = get_synodic_period(ipl);
2589 tjdcon = tjd0 + ((floor) ((tjd_start - tjd0) / dsynperiod) + 1) * dsynperiod;
2590 ds = 100;
2591 while (ds > 0.5) {
2592 if (swe_calc(tjdcon, ipl, epheflag|SEFLG_SPEED, x, serr) == ERR)
2593 return ERR;
2594 if (swe_calc(tjdcon, SE_SUN, epheflag|SEFLG_SPEED, xs, serr) == ERR)
2595 return ERR;
2596 ds = swe_degnorm(x[0] - xs[0] - daspect);
2597 if (ds > 180) ds -= 360;
2598 tjdcon -= ds / (x[3] - xs[3]);
2599 }
2600 *tjd = tjdcon;
2601 return OK;
2602 }
2603
get_asc_obl_with_sun(double tjd_start,int32 ipl,char * star,int32 helflag,int32 evtyp,double dperiod,double * dgeo,double * tjdret,char * serr)2604 static int32 get_asc_obl_with_sun(double tjd_start, int32 ipl, char *star, int32 helflag, int32 evtyp, double dperiod, double *dgeo, double *tjdret, char *serr)
2605 {
2606 int i, retval;
2607 int32 is_acronychal = FALSE;
2608 int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
2609 double dsunpl = 1, dsunpl_save, dsunpl_test, tjd, daystep;
2610 AS_BOOL desc_obl = FALSE, retro = FALSE;
2611 if (evtyp == SE_EVENING_LAST || evtyp == SE_EVENING_FIRST)
2612 desc_obl = TRUE;
2613 if (evtyp == SE_MORNING_FIRST || evtyp == SE_EVENING_LAST)
2614 retro = TRUE;
2615 if (evtyp == SE_ACRONYCHAL_RISING)
2616 desc_obl = TRUE;
2617 if (evtyp == SE_ACRONYCHAL_RISING || evtyp == SE_ACRONYCHAL_SETTING) {
2618 is_acronychal = TRUE;
2619 if (ipl != SE_MOON)
2620 retro = TRUE;
2621 }
2622 // if (evtyp == 3 || evtyp == 4)
2623 // dangsearch = 180;
2624 /* find date when sun and object have the same ascensio obliqua */
2625 tjd = tjd_start;
2626 dsunpl_save = -999999999;
2627 retval = get_asc_obl_diff(tjd, ipl, star, epheflag, dgeo, desc_obl, is_acronychal, &dsunpl, serr);
2628 if (retval != OK) /* retval may be ERR or -2 */
2629 return retval;
2630 daystep = 20;
2631 i = 0;
2632 while (dsunpl_save == -999999999 ||
2633 /*fabs(dsunpl - dsunpl_save) > 180 ||*/
2634 fabs(dsunpl) + fabs(dsunpl_save) > 180 ||
2635 (retro && !(dsunpl_save < 0 && dsunpl >= 0)) ||
2636 (!retro && !(dsunpl_save >= 0 && dsunpl < 0))) {
2637 i++;
2638 if (i > 5000) {
2639 sprintf(serr, "loop in get_asc_obl_with_sun() (1)");
2640 return ERR;
2641 }
2642 dsunpl_save = dsunpl;
2643 tjd += 10.0;
2644 if (dperiod > 0 && tjd - tjd_start > dperiod)
2645 return -2;
2646 retval = get_asc_obl_diff(tjd, ipl, star, epheflag, dgeo, desc_obl, is_acronychal, &dsunpl, serr);
2647 if (retval != OK) /* retval may be ERR or -2 */
2648 return retval;
2649 }
2650 tjd_start = tjd - daystep;
2651 daystep /= 2.0;
2652 tjd = tjd_start + daystep;
2653 retval = get_asc_obl_diff(tjd, ipl, star, epheflag, dgeo, desc_obl, is_acronychal, &dsunpl_test, serr);
2654 if (retval != OK) /* retval may be ERR or -2 */
2655 return retval;
2656 i = 0;
2657 while (fabs(dsunpl) > 0.00001) {
2658 i++;
2659 if (i > 5000) {
2660 sprintf(serr, "loop in get_asc_obl_with_sun() (2)");
2661 return ERR;
2662 }
2663 if (dsunpl_save * dsunpl_test >= 0) {
2664 dsunpl_save = dsunpl_test;
2665 tjd_start = tjd;
2666 } else {
2667 dsunpl = dsunpl_test;
2668 }
2669 daystep /= 2.0;
2670 tjd = tjd_start + daystep;
2671 retval = get_asc_obl_diff(tjd, ipl, star, epheflag, dgeo, desc_obl, is_acronychal, &dsunpl_test, serr);
2672 if (retval != OK) /* retval may be ERR or -2 */
2673 return retval;
2674 }
2675 *tjdret = tjd;
2676 return OK;
2677 }
2678
2679 #if 0
2680 /* works only for fixed stars */
2681 static int32 get_asc_obl_with_sun_old(double tjd_start, int32 ipl, char *star, int32 helflag, int32 TypeEvent, double *dgeo, double *tjdret, char *serr)
2682 {
2683 int retval;
2684 int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
2685 double dsunpl = 1, tjd, daystep, dsunpl_save;
2686 double dsynperiod = 367;
2687 double dangsearch = 0;
2688 AS_BOOL desc_obl = FALSE;
2689 if (TypeEvent == 2 || TypeEvent == 3)
2690 desc_obl = TRUE;
2691 if (TypeEvent == 3 || TypeEvent == 4)
2692 dangsearch = 180;
2693 /* find date when sun and object have the same ascensio obliqua */
2694 daystep = dsynperiod;
2695 tjd = tjd_start;
2696 retval = get_asc_obl_diff(tjd, ipl, star, epheflag, dgeo, desc_obl, &dsunpl, serr);
2697 if (retval != OK) /* retval may be ERR or -2 */
2698 return retval;
2699 while (dsunpl < 359.99999) {
2700 dsunpl_save = dsunpl;
2701 daystep /= 2.0;
2702 retval = get_asc_obl_diff(tjd + daystep, ipl, star, epheflag, dgeo, desc_obl, &dsunpl, serr);
2703 if (retval != OK) /* retval may be ERR or -2 */
2704 return retval;
2705 if (dsunpl > dsunpl_save)
2706 tjd += daystep;
2707 else
2708 dsunpl = dsunpl_save;
2709 }
2710 *tjdret = tjd;
2711 return OK;
2712 }
2713 #endif
2714
2715 #if 0
2716 /* works only for fixed stars */
2717 static int32 get_asc_obl_acronychal(double tjd_start, int32 ipl, char *star, int32 helflag, int32 TypeEvent, double *dgeo, double *tjdret, char *serr)
2718 {
2719 int retval;
2720 int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
2721 double dsunpl = 1, tjd, daystep, dsunpl_save;
2722 double dsynperiod = 367;
2723 double aosun, aopl;
2724 AS_BOOL sun_desc = TRUE, obj_desc = FALSE;
2725 daystep = dsynperiod;
2726 tjd = tjd_start;
2727 if (TypeEvent == 4) {
2728 sun_desc = FALSE;
2729 obj_desc = TRUE;
2730 }
2731 /* ascensio (descensio) obliqua of sun */
2732 retval = get_asc_obl(tjd, SE_SUN, "", epheflag, dgeo, sun_desc, &aosun, serr);
2733 if (retval != OK) /* retval may be ERR or -2 */
2734 return retval;
2735 /* ascensio (descensio) obliqua of body */
2736 retval = get_asc_obl(tjd, ipl, star, epheflag, dgeo, obj_desc, &aopl, serr);
2737 if (retval != OK) /* retval may be ERR or -2 */
2738 return retval;
2739 dsunpl = swe_degnorm(aosun - aopl + 180);
2740 while (dsunpl < 359.99999) {
2741 dsunpl_save = dsunpl;
2742 daystep /= 2.0;
2743 /* ascensio (descensio) obliqua of sun */
2744 retval = get_asc_obl(tjd+daystep, SE_SUN, "", epheflag, dgeo, sun_desc, &aosun, serr);
2745 if (retval != OK) /* retval may be ERR or -2 */
2746 return retval;
2747 /* ascensio (descensio) obliqua of body */
2748 retval = get_asc_obl(tjd+daystep, ipl, star, epheflag, dgeo, obj_desc, &aopl, serr);
2749 if (retval != OK) /* retval may be ERR or -2 */
2750 return retval;
2751 dsunpl = swe_degnorm(aosun - aopl + 180);
2752 if (dsunpl > dsunpl_save)
2753 tjd += daystep;
2754 else
2755 dsunpl = dsunpl_save;
2756 }
2757 *tjdret = tjd;
2758 return OK;
2759 }
2760 #endif
2761
get_heliacal_day(double tjd,double * dgeo,double * datm,double * dobs,char * ObjectName,int32 helflag,int32 TypeEvent,double * thel,char * serr)2762 static int32 get_heliacal_day(double tjd, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 helflag, int32 TypeEvent, double *thel, char *serr)
2763 {
2764 int32 i, visible_at_sunsetrise, is_rise_or_set = 0, ndays, retval, retval_old;
2765 double direct_day = 0, direct_time = 0, tfac, tend, daystep, tday, vdelta, tret;
2766 double darr[30], vd, dmag, div;
2767 int32 ipl = DeterObject(ObjectName);
2768 /*
2769 * find the day and minute on which the object becomes visible
2770 */
2771 switch (TypeEvent) {
2772 /* morning first */
2773 case 1: is_rise_or_set = SE_CALC_RISE;
2774 direct_day = 1; direct_time = -1;
2775 break;
2776 /* evening last */
2777 case 2: is_rise_or_set = SE_CALC_SET;
2778 direct_day = -1; direct_time = 1;
2779 break;
2780 /* evening first */
2781 case 3: is_rise_or_set = SE_CALC_SET;
2782 direct_day = 1; direct_time = 1;
2783 break;
2784 /* morning last */
2785 case 4: is_rise_or_set = SE_CALC_RISE;
2786 direct_day = -1; direct_time = -1;
2787 break;
2788 }
2789 tfac = 1;
2790 switch (ipl) {
2791 case SE_MOON:
2792 ndays = 16;
2793 daystep = 1;
2794 break;
2795 case SE_MERCURY:
2796 ndays = 60; tjd -= 0 * direct_day;
2797 daystep = 5;
2798 tfac = 5;
2799 break;
2800 case SE_VENUS:
2801 ndays = 300; tjd -= 30 * direct_day;
2802 daystep = 5;
2803 if (TypeEvent >= 3) {
2804 daystep = 15;
2805 tfac = 3;
2806 }
2807 break;
2808 case SE_MARS:
2809 ndays = 400;
2810 daystep = 15;
2811 tfac = 5;
2812 break;
2813 case SE_SATURN:
2814 ndays = 300;
2815 daystep = 20;
2816 tfac = 5;
2817 break;
2818 case -1:
2819 ndays = 300;
2820 if (call_swe_fixstar_mag(ObjectName, &dmag, serr) == ERR) {
2821 return ERR;
2822 }
2823 daystep = 15;
2824 tfac = 10;
2825 if (dmag > 2) {
2826 daystep = 15;
2827 }
2828 if (dmag < 0) {
2829 tfac = 3;
2830 }
2831 break;
2832 default:
2833 ndays = 300;
2834 daystep = 15;
2835 tfac = 3;
2836 break;
2837 }
2838 tend = tjd + ndays * direct_day;
2839 retval_old = -2;
2840 for (tday = tjd, i = 0;
2841 (direct_day > 0 && tday < tend) || (direct_day < 0 && tday > tend);
2842 tday += daystep * direct_day, i++) {
2843 vdelta = -100;
2844 if (i > 0)
2845 tday -= 0.3 * direct_day;
2846 if ((retval = my_rise_trans(tday, SE_SUN, "", is_rise_or_set, helflag, dgeo, datm, &tret, serr)) == ERR) {
2847 return ERR;
2848 }
2849 /* sun does not rise: try next day */
2850 if (retval == -2) {
2851 retval_old = retval;
2852 continue;
2853 }
2854 retval = swe_vis_limit_mag(tret, dgeo, datm, dobs, ObjectName, helflag, darr, serr);
2855 if (retval == ERR)
2856 return ERR;
2857 #if 1
2858 /* object has appeared above horizon: reduce daystep */
2859 if (retval_old == -2 && retval >= 0 && daystep > 1) {
2860 retval_old = retval;
2861 tday -= daystep * direct_day;
2862 daystep = 1;
2863 /* Note: beyond latitude 55N (?), Mars can have a morning last.
2864 * If the period of visibility is less than 5 days, we may miss the
2865 * event. I don't know if this happens */
2866 if (ipl >= SE_MARS || ipl == -1)
2867 daystep = 5;
2868 continue;
2869 }
2870 retval_old = retval;
2871 #endif
2872 /* object below horizon: try next day */
2873 if (retval == -2)
2874 continue;
2875 vdelta = darr[0] - darr[7];
2876 /* find minute of object's becoming visible */
2877 div = 1440.0;
2878 // div = 86400.0;
2879 // printf("tret=%f, m1=%f, m2=%f\n", tret, darr[0], darr[1]);
2880 vd = -1;
2881 visible_at_sunsetrise = 1;
2882 while (retval != -2 && (vd = darr[0] - darr[7]) < 0) {
2883 visible_at_sunsetrise = 0;
2884 if (vd < -1.0)
2885 tret += 5.0 / div * direct_time * tfac;
2886 else if (vd < -0.5)
2887 tret += 2.0 / div * direct_time * tfac;
2888 else if (vd < -0.1)
2889 tret += 1.0 / div * direct_time * tfac;
2890 else
2891 tret += 1.0 / div * direct_time;
2892 retval = swe_vis_limit_mag(tret, dgeo, datm, dobs, ObjectName, helflag, darr, serr);
2893 if (retval == ERR)
2894 return ERR;
2895 }
2896 /* if possible move a bit away from sunset, where vis_limit_mag() has strange behaviour */
2897 if (visible_at_sunsetrise) {
2898 for (i = 0; i < 10; i++) {
2899 if ((retval = swe_vis_limit_mag(tret + 1.0 / div * direct_time, dgeo, datm, dobs, ObjectName, helflag, darr, serr)) >= 0
2900 && darr[0] - darr[7] > vd) {
2901 vd = darr[0] - darr[7];
2902 tret += 1.0 / div * direct_time;
2903 }
2904 }
2905 }
2906 vdelta = darr[0] - darr[7];
2907 /* object is visible, save time of appearance */
2908 if (vdelta > 0) {
2909 if ((ipl >= SE_MARS || ipl == -1) && daystep > 1) {
2910 tday -= daystep * direct_day;
2911 daystep = 1;
2912 } else {
2913 *thel = tret;
2914 return OK;
2915 }
2916 }
2917 }
2918 sprintf(serr, "heliacal event does not happen");
2919 //printf("%s\n", serr);
2920 return -2;
2921 }
2922
time_optimum_visibility(double tjd,double * dgeo,double * datm,double * dobs,char * ObjectName,int32 helflag,double * tret,char * serr)2923 static int32 time_optimum_visibility(double tjd, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 helflag, double *tret, char *serr)
2924 {
2925 int32 retval, retval_sv, i;
2926 double t1, t2, vl1, vl2, d, darr[10], phot_scot_opic, phot_scot_opic_sv;
2927 // double vl;
2928 int t_has_changed;
2929 *tret = tjd;
2930 retval = swe_vis_limit_mag(tjd, dgeo, datm, dobs, ObjectName, helflag, darr, serr);
2931 if (retval == ERR) return ERR;
2932 retval_sv = retval;
2933 //vl = darr[0] - darr[7];
2934 //vl = -1;
2935 t1 = tjd;
2936 t2 = tjd;
2937 vl1 = -1;
2938 vl2 = -1;
2939 //printf("begin tret=%f, dvl=%f\n", tjd, darr[0] - darr[7]);
2940 phot_scot_opic_sv = retval & SE_SCOTOPIC_FLAG;
2941 for (i = 0, d = 100.0 / 86400.0; i < 3; i++, d /= 10.0) {
2942 // fprintf(stderr, "i= %d\n", i);
2943 t1 += d;
2944 t_has_changed = 0;
2945 while((retval = swe_vis_limit_mag(t1 - d, dgeo, datm, dobs, ObjectName, helflag, darr, serr)) >= 0
2946 && darr[0] > darr[7]
2947 && darr[0] - darr[7] > vl1) {
2948 t1 -= d; vl1 = darr[0] - darr[7];
2949 t_has_changed = 1;
2950 //fprintf(stderr, "vl1=%f %d vlm=%f, obm=%f, t=%f\n", vl, retval, darr[0], darr[7], tjd + d);
2951 retval_sv = retval;
2952 phot_scot_opic_sv = retval & SE_SCOTOPIC_FLAG;
2953 /* printf("1: %f\n", darr[8]);*/
2954 }
2955 if (t_has_changed == 0)
2956 t1 -= d; /* revert initial addition */
2957 if (retval == ERR) return ERR;
2958 }
2959 for (i = 0, d = 100.0 / 86400.0; i < 3; i++, d /= 10.0) {
2960 t2 -= d;
2961 t_has_changed = 0;
2962 while((retval = swe_vis_limit_mag(t2 + d, dgeo, datm, dobs, ObjectName, helflag, darr, serr)) >= 0
2963 && darr[0] > darr[7]
2964 && darr[0] - darr[7] > vl2) {
2965 t2 += d; vl2 = darr[0] - darr[7];
2966 t_has_changed = 1;
2967 //fprintf(stderr, "vl2=%f %d vlm=%f, obm=%f, t=%f\n", vl, retval, darr[0], darr[7], tjd + d);
2968 retval_sv = retval;
2969 phot_scot_opic_sv = retval & SE_SCOTOPIC_FLAG;
2970 /* printf("2: %f\n", darr[8]);*/
2971 }
2972 if (t_has_changed == 0)
2973 t2 += d; /* revert initial subtraction */
2974 if (retval == ERR) return ERR;
2975 }
2976 if (vl2 > vl1)
2977 tjd = t2;
2978 else
2979 tjd = t1;
2980 /* printf("3: %f <-> %f\n", darr[8], phot_scot_opic_sv);*/
2981 *tret = tjd;
2982 //printf("end tret=%f, dvl=%f\n", tjd, darr[0] - darr[7]);
2983 if (retval >= 0) {
2984 /* search for optimum came to an end because change scotopic/photopic: */
2985 phot_scot_opic = (retval & SE_SCOTOPIC_FLAG);
2986 if (phot_scot_opic_sv != phot_scot_opic) {
2987 /* calling function writes warning into serr */
2988 printf ("hallo -2\n");
2989 return -2;
2990 }
2991 /* valid result found but it is close to the scotopic/photopic limit */
2992 if (retval_sv & SE_MIXEDOPIC_FLAG) {
2993 printf ("hallo -2\n");
2994 return -2;
2995 }
2996 }
2997 return OK;
2998 }
2999
time_limit_invisible(double tjd,double * dgeo,double * datm,double * dobs,char * ObjectName,int32 helflag,int32 direct,double * tret,char * serr)3000 static int32 time_limit_invisible(double tjd, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 helflag, int32 direct, double *tret, char *serr)
3001 {
3002 int32 retval, retval_sv, i, ncnt = 3;
3003 double d = 0, darr[10], phot_scot_opic, phot_scot_opic_sv;
3004 double d0 = 100.0 / 86400.0;
3005 *tret = tjd;
3006 if (strcmp(ObjectName, "moon") == 0) {
3007 d0 *= 10;
3008 ncnt = 4;
3009 }
3010 retval = swe_vis_limit_mag(tjd + d * direct, dgeo, datm, dobs, ObjectName, helflag, darr, serr);
3011 if (retval == ERR) return ERR;
3012 retval_sv = retval;
3013 phot_scot_opic_sv = retval & SE_SCOTOPIC_FLAG;
3014 for (i = 0, d = d0; i < ncnt; i++, d /= 10.0) {
3015 while((retval = swe_vis_limit_mag(tjd + d * direct, dgeo, datm, dobs, ObjectName, helflag, darr, serr)) >= 0
3016 && darr[0] > darr[7]) {
3017 tjd += d * direct;
3018 retval_sv = retval;
3019 phot_scot_opic_sv = retval & SE_SCOTOPIC_FLAG;
3020 /* printf("%d: %f\n", direct, darr[8]); */
3021 }
3022 }
3023 /* printf("4: %f, %f/%f %f <-> %f\n", darr[8], darr[0], darr[7], tjd, phot_scot_opic_sv); */
3024 *tret = tjd;
3025 /* if object disappears at setting, retval is -2, but we want it OK, and
3026 * also suppress the warning "object is below local horizon" */
3027 *serr = '\0';
3028 if (retval >= 0) {
3029 /* search for limit came to an end because change scotopic/photopic: */
3030 phot_scot_opic = (retval & SE_SCOTOPIC_FLAG);
3031 if (phot_scot_opic_sv != phot_scot_opic) {
3032 /* calling function writes warning into serr */
3033 return -2;
3034 }
3035 /* valid result found but it is close to the scotopic/photopic limit */
3036 if (retval_sv & SE_MIXEDOPIC_FLAG) {
3037 return -2;
3038 }
3039 }
3040 return OK;
3041 }
3042
get_acronychal_day(double tjd,double * dgeo,double * datm,double * dobs,char * ObjectName,int32 helflag,int32 TypeEvent,double * thel,char * serr)3043 static int32 get_acronychal_day(double tjd, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 helflag, int32 TypeEvent, double *thel, char *serr) {
3044 double tret, tret_dark, darr[30], dtret;
3045 /* x[6], xaz[6], alto, azio, alto_dark, azio_dark;*/
3046 int32 retval, is_rise_or_set, direct;
3047 int32 ipl = DeterObject(ObjectName);
3048 helflag |= SE_HELFLAG_VISLIM_PHOTOPIC;
3049 /*int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);*/
3050 /* int32 iflag = epheflag | SEFLG_EQUATORIAL | SEFLG_TOPOCTR;*/
3051 if (TypeEvent == 3 || TypeEvent == 5) {
3052 is_rise_or_set = SE_CALC_RISE;
3053 /* tret = tjdc - 3;
3054 if (ipl >= SE_MARS)
3055 tret = tjdc - 3;*/
3056 direct = -1;
3057 } else {
3058 is_rise_or_set = SE_CALC_SET;
3059 /*tret = tjdc + 3;
3060 if (ipl >= SE_MARS)
3061 tret = tjdc + 3;*/
3062 direct = 1;
3063 }
3064 dtret = 999;
3065 #if 0
3066 while (fabs(dtret) > 0.5) {
3067 #else
3068 while (fabs(dtret) > 0.5 / 1440.0) {
3069 #endif
3070 tjd += 0.7 * direct;
3071 if (direct < 0) tjd -= 1;
3072 retval = my_rise_trans(tjd, ipl, ObjectName, is_rise_or_set, helflag, dgeo, datm, &tjd, serr);
3073 if (retval == ERR) return ERR;
3074 retval = swe_vis_limit_mag(tjd, dgeo, datm, dobs, ObjectName, helflag, darr, serr);
3075 if (retval == ERR) return ERR;
3076 while(darr[0] < darr[7]) {
3077 tjd += 10.0 / 1440.0 * -direct;
3078 retval = swe_vis_limit_mag(tjd, dgeo, datm, dobs, ObjectName, helflag, darr, serr);
3079 if (retval == ERR) return ERR;
3080 }
3081 retval = time_limit_invisible(tjd, dgeo, datm, dobs, ObjectName, helflag | SE_HELFLAG_VISLIM_DARK, direct, &tret_dark, serr);
3082 if (retval == ERR) return ERR;
3083 retval = time_limit_invisible(tjd, dgeo, datm, dobs, ObjectName, helflag | SE_HELFLAG_VISLIM_NOMOON, direct, &tret, serr);
3084 if (retval == ERR) return ERR;
3085 #if 0
3086 if (azalt_cart(tret_dark, dgeo, datm, ObjectName, helflag, darr, serr) == ERR)
3087 return ERR;
3088 if (azalt_cart(tret, dgeo, datm, ObjectName, helflag, darr+6, serr) == ERR)
3089 return ERR;
3090 dtret = acos(swi_dot_prod_unit(darr+3, darr+9)) / DEGTORAD;
3091 #else
3092 dtret = fabs(tret - tret_dark);
3093 #endif
3094 }
3095 if (azalt_cart(tret, dgeo, datm, "sun", helflag, darr, serr) == ERR)
3096 return ERR;
3097 *thel = tret;
3098 if (darr[1] < -12) {
3099 sprintf(serr, "acronychal rising/setting not available, %f", darr[1]);
3100 return OK;
3101 } else {
3102 sprintf(serr, "solar altitude, %f", darr[1]);
3103 }
3104 return OK;
3105 }
3106
3107 static int32 get_heliacal_details(double tday, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 TypeEvent, int32 helflag, double *dret, char *serr)
3108 {
3109 int32 i, retval, direct;
3110 AS_BOOL optimum_undefined, limit_1_undefined, limit_2_undefined;
3111 /* find next optimum visibility */
3112 optimum_undefined = FALSE;
3113 retval = time_optimum_visibility(tday, dgeo, datm, dobs, ObjectName, helflag, &(dret[1]), serr);
3114 if (retval == ERR) return ERR;
3115 if (retval == -2) {
3116 retval = OK;
3117 optimum_undefined = TRUE; /* change photopic <-> scotopic vision */
3118 }
3119 /* find moment of becoming visible */
3120 direct = 1;
3121 if (TypeEvent == 1 || TypeEvent == 4)
3122 direct = -1;
3123 limit_1_undefined = FALSE;
3124 retval = time_limit_invisible(tday, dgeo, datm, dobs, ObjectName, helflag, direct, &(dret[0]), serr);
3125 if (retval == ERR) return ERR;
3126 if (retval == -2) {
3127 retval = OK;
3128 limit_1_undefined = TRUE; /* change photopic <-> scotopic vision */
3129 }
3130 /* find moment of end of visibility */
3131 direct *= -1;
3132 limit_2_undefined = FALSE;
3133 retval = time_limit_invisible(dret[1], dgeo, datm, dobs, ObjectName, helflag, direct, &(dret[2]), serr);
3134 if (retval == ERR) return ERR;
3135 if (retval == -2) {
3136 retval = OK;
3137 limit_2_undefined = TRUE; /* change photopic <-> scotopic vision */
3138 }
3139 /* correct sequence of times:
3140 * with event types 2 and 3 swap dret[0] and dret[2] */
3141 if (TypeEvent == 2 || TypeEvent == 3) {
3142 tday = dret[2];
3143 dret[2] = dret[0];
3144 dret[0] = tday;
3145 i = (int) limit_1_undefined;
3146 limit_1_undefined = limit_2_undefined;
3147 limit_2_undefined = (AS_BOOL) i;
3148 }
3149 /*if (retval == OK && dret[0] == dret[1]) */
3150 if (optimum_undefined || limit_1_undefined || limit_2_undefined) {
3151 sprintf(serr, "return values [");
3152 if (limit_1_undefined)
3153 strcat(serr, "0,");
3154 if (optimum_undefined)
3155 strcat(serr, "1,");
3156 if (limit_2_undefined)
3157 strcat(serr, "2,");
3158 strcat(serr, "] are uncertain due to change between photopic and scotopic vision");
3159 }
3160 return OK;
3161 }
3162
3163 static int32 heliacal_ut_vis_lim(double tjd_start, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 TypeEventIn, int32 helflag, double *dret, char *serr_ret)
3164 {
3165 int i;
3166 double d, darr[10], direct = 1, tjd, tday;
3167 int32 epheflag, retval = OK, helflag2;
3168 int32 iflag, ipl;
3169 int32 TypeEvent = TypeEventIn;
3170 char serr[AS_MAXCH];
3171 for (i = 0; i < 10; i++)
3172 dret[i] = 0;
3173 *dret = tjd_start;
3174 *serr = '\0';
3175 ipl = DeterObject(ObjectName);
3176 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
3177 iflag = SEFLG_TOPOCTR | SEFLG_EQUATORIAL | epheflag;
3178 if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
3179 iflag |= SEFLG_NONUT | SEFLG_TRUEPOS;
3180 if (ipl == SE_MERCURY)
3181 tjd = tjd_start - 30;
3182 else
3183 tjd = tjd_start - 50; /* -50 makes sure, that no event is missed,
3184 * but may return an event before start date */
3185 helflag2 = helflag;
3186 /*helflag2 &= ~SE_HELFLAG_HIGH_PRECISION;*/
3187 /*
3188 * heliacal event
3189 */
3190 if (ipl == SE_MERCURY || ipl == SE_VENUS || TypeEvent <= 2) {
3191 if (ipl == -1) {
3192 /* find date when star rises with sun (cosmic rising) */
3193 retval = get_asc_obl_with_sun(tjd, ipl, ObjectName, helflag, TypeEvent, 0, dgeo, &tjd, serr);
3194 if (retval != OK)
3195 goto swe_heliacal_err; /* retval may be -2 or ERR */
3196 } else {
3197 /* find date of conjunction of object with sun */
3198 if ((retval = find_conjunct_sun(tjd, ipl, helflag, TypeEvent, &tjd, serr)) == ERR) {
3199 goto swe_heliacal_err;
3200 }
3201 }
3202 /* find the day and minute on which the object becomes visible */
3203 retval = get_heliacal_day(tjd, dgeo, datm, dobs, ObjectName, helflag2, TypeEvent, &tday, serr);
3204 if (retval != OK)
3205 goto swe_heliacal_err;
3206 /*
3207 * acronychal event
3208 */
3209 } else {
3210 if (1 || ipl == -1) {
3211 /*retval = get_asc_obl_acronychal(tjd, ipl, ObjectName, helflag2, TypeEvent, dgeo, &tjd, serr);*/
3212 retval = get_asc_obl_with_sun(tjd, ipl, ObjectName, helflag, TypeEvent, 0, dgeo, &tjd, serr);
3213 if (retval != OK)
3214 goto swe_heliacal_err;
3215 } else {
3216 /* find date of conjunction of object with sun */
3217 if ((retval = find_conjunct_sun(tjd, ipl, helflag, TypeEvent, &tjd, serr)) == ERR)
3218 goto swe_heliacal_err;
3219 }
3220 tday = tjd;
3221 retval = get_acronychal_day(tjd, dgeo, datm, dobs, ObjectName, helflag2, TypeEvent, &tday, serr);
3222 if (retval != OK)
3223 goto swe_heliacal_err;
3224 }
3225 dret[0] = tday;
3226 if (!(helflag & SE_HELFLAG_NO_DETAILS)) {
3227 /* more precise event times for
3228 * - morning first, evening last
3229 * - venus and mercury's evening first and morning last
3230 */
3231 if (ipl == SE_MERCURY || ipl == SE_VENUS || TypeEvent <= 2) {
3232 retval = get_heliacal_details(tday, dgeo, datm, dobs, ObjectName, TypeEvent, helflag2, dret, serr);
3233 if (retval == ERR) goto swe_heliacal_err;
3234 } else if ((0)) {
3235 if (TypeEvent == 4 || TypeEvent == 6) direct = -1;
3236 for (i = 0, d = 100.0 / 86400.0; i < 3; i++, d /= 10.0) {
3237 while((retval = swe_vis_limit_mag(*dret + d * direct, dgeo, datm, dobs, ObjectName, helflag, darr, serr)) == -2 || (retval >= 0 && darr[0] < darr[7])) {
3238 *dret += d * direct;
3239 }
3240 }
3241 /* the last time step must be added */
3242 if (retval == OK)
3243 *dret += 1.0 / 86400.0 * direct;
3244 }
3245 } /* if (1) */
3246 swe_heliacal_err:
3247 if (serr_ret != NULL && *serr != '\0')
3248 strcpy(serr_ret, serr);
3249 return retval;
3250 }
3251
3252 /*###################################################################*/
3253 static int32 moon_event_vis_lim(double tjdstart, double *dgeo, double *datm, double *dobs, int32 TypeEvent, int32 helflag, double *dret, char *serr_ret)
3254 {
3255 double tjd, trise;
3256 char serr[AS_MAXCH];
3257 char ObjectName[30];
3258 int32 iflag, ipl, retval, helflag2, direct;
3259 int32 epheflag = helflag & (SEFLG_JPLEPH|SEFLG_SWIEPH|SEFLG_MOSEPH);
3260 dret[0] = tjdstart; /* will be returned in error case */
3261 if (TypeEvent == 1 || TypeEvent == 2) {
3262 if (serr_ret != NULL)
3263 strcpy(serr_ret, "error: the moon has no morning first or evening last");
3264 return ERR;
3265 }
3266 strcpy(ObjectName, "moon");
3267 ipl = SE_MOON;
3268 iflag = SEFLG_TOPOCTR | SEFLG_EQUATORIAL | epheflag;
3269 if (!(helflag & SE_HELFLAG_HIGH_PRECISION))
3270 iflag |= SEFLG_NONUT|SEFLG_TRUEPOS;
3271 helflag2 = helflag;
3272 helflag2 &= ~SE_HELFLAG_HIGH_PRECISION;
3273 /* check Synodic/phase Period */
3274 tjd = tjdstart - 30; /* -50 makes sure, that no event is missed,
3275 * but may return an event before start date */
3276 if ((retval = find_conjunct_sun(tjd, ipl, helflag, TypeEvent, &tjd, serr)) == ERR)
3277 return ERR;
3278 /* find the day and minute on which the object becomes visible */
3279 retval = get_heliacal_day(tjd, dgeo, datm, dobs, ObjectName, helflag2, TypeEvent, &tjd, serr);
3280 if (retval != OK)
3281 goto moon_event_err;
3282 dret[0] = tjd;
3283 /* find next optimum visibility */
3284 retval = time_optimum_visibility(tjd, dgeo, datm, dobs, ObjectName, helflag, &tjd, serr);
3285 if (retval == ERR) goto moon_event_err;
3286 dret[1] = tjd;
3287 /* find moment of becoming visible */
3288 /* Note: On the day of first light the moon may become visible
3289 * already during day. It also may appear during day, disappear again
3290 * and then reappear after sunset */
3291 direct = 1;
3292 if (TypeEvent == 4)
3293 direct = -1;
3294 retval = time_limit_invisible(tjd, dgeo, datm, dobs, ObjectName, helflag, direct, &tjd, serr);
3295 if (retval == ERR) goto moon_event_err;
3296 dret[2] = tjd;
3297 /* find moment of end of visibility */
3298 direct *= -1;
3299 retval = time_limit_invisible(dret[1], dgeo, datm, dobs, ObjectName, helflag, direct, &tjd, serr);
3300 dret[0] = tjd;
3301 if (retval == ERR) goto moon_event_err;
3302 #if 1
3303 /* if the moon is visible before sunset, we return sunset as start time */
3304 if (TypeEvent == 3) {
3305 if ((retval = my_rise_trans(tjd, SE_SUN, "", SE_CALC_SET, helflag, dgeo, datm, &trise, serr)) == ERR)
3306 return ERR;
3307 if (trise < dret[1]) {
3308 dret[0] = trise;
3309 /* do not warn, it happens too often */
3310 /*strcpy(serr, "start time given is sunset, but moon is observable before that");*/
3311 }
3312 /* if the moon is visible after sunrise, we return sunrise as end time */
3313 } else {
3314 if ((retval = my_rise_trans(dret[1], SE_SUN, "", SE_CALC_RISE, helflag, dgeo, datm, &trise, serr)) == ERR)
3315 return ERR;
3316 if (dret[0] > trise) {
3317 dret[0] = trise;
3318 /* do not warn, it happens too often */
3319 /*strcpy(serr, "end time given is sunrise, but moon is observable after that");*/
3320 }
3321 }
3322 #endif
3323 /* correct order of the three times: */
3324 if (TypeEvent == 4) {
3325 tjd = dret[0];
3326 dret[0] = dret[2];
3327 dret[2] = tjd;
3328 }
3329 moon_event_err:
3330 if (serr_ret != NULL && *serr != '\0')
3331 strcpy(serr_ret, serr);
3332 return retval;
3333 }
3334
3335 static int32 MoonEventJDut(double JDNDaysUTStart, double *dgeo, double *datm, double *dobs, int32 TypeEvent, int32 helflag, double *dret, char *serr)
3336 {
3337 int32 avkind = helflag & SE_HELFLAG_AVKIND;
3338 if (avkind)
3339 return moon_event_arc_vis(JDNDaysUTStart, dgeo, datm, dobs, TypeEvent, helflag, dret, serr);
3340 else
3341 return moon_event_vis_lim(JDNDaysUTStart, dgeo, datm, dobs, TypeEvent, helflag, dret, serr);
3342 }
3343
3344 static int32 heliacal_ut(double JDNDaysUTStart, double *dgeo, double *datm, double *dobs, char *ObjectName, int32 TypeEventIn, int32 helflag, double *dret, char *serr_ret)
3345 {
3346 int32 avkind = helflag & SE_HELFLAG_AVKIND;
3347 if (avkind)
3348 return heliacal_ut_arc_vis(JDNDaysUTStart, dgeo, datm, dobs, ObjectName, TypeEventIn, helflag, dret, serr_ret);
3349 else
3350 return heliacal_ut_vis_lim(JDNDaysUTStart, dgeo, datm, dobs, ObjectName, TypeEventIn, helflag, dret, serr_ret);
3351 }
3352
3353 /*' Magn [-]
3354 ' tjd_ut start date (JD) for event search
3355 ' dgeo[3] geogr. longitude, latitude, eye height (m above sea level)
3356 ' datm[4] atm. pressure, temperature, RH, and VR
3357 ' - pressure atmospheric pressure (mbar, =hPa) default 1013.25hPa
3358 ' - temperature deg C, default 15 deg C (if at
3359 ' If both attemp and atpress are 0, a temperature and
3360 ' atmospheric pressure are estimated from the above-mentioned
3361 ' default values and the height above sea level.
3362 ' - RH relative humidity in %
3363 ' - VR VR>=1: the Meteorological range: default 40 km
3364 ' 1>VR>0: the ktot (so the total atmospheric coefficient):
3365 ' a good default would be 0.25
3366 ' VR=-1: the ktot is calculated from the other atmospheric
3367 ' constants.
3368 '
3369 ' dobs[6] observer parameters
3370 ' - age [Year] default 36, experienced sky observer in ancient times
3371 ' optimum age is 23
3372 ' - SN Snellen factor of the visual aquity of the observer
3373 ' default 1
3374 ' see: http://www.i-see.org/eyecharts.html#make-your-own
3375 ' The following parameters of dobs[] are only relevant if the flag
3376 ' SE_HELFLAG_OPTICAL_PARAMS is set:
3377 ' - is_binocular 0 = monocular, 1 = binocular (actually a boolean)
3378 ' - OpticMagn telescope magnification:
3379 ' 0 = default to naked eye (binocular), 1 = naked eye
3380 ' - OpticDia optical aperture (telescope diameter) in mm
3381 ' - OpticTrans optical transmission
3382 '
3383 ' TypeEvent 1 morning first
3384 ' 2 evening last
3385 ' 3 evening first
3386 ' 4 morning last
3387 ' dret output: time (tjd_ut) of heliacal event
3388 ' dret[0]: beginning of visibility (Julian day number)
3389 ' dret[1]: optimum visibility (Julian day number; 0 if SE_HELFLAG_AV)
3390 ' dret[2]: end of visibility (Julian day number; 0 if SE_HELFLAG_AV)
3391 ' see http://www.iol.ie/~geniet/eng/atmoastroextinction.htm
3392 */
3393 int32 CALL_CONV swe_heliacal_ut(double JDNDaysUTStart, double *dgeo, double *datm, double *dobs, char *ObjectNameIn, int32 TypeEvent, int32 helflag, double *dret, char *serr_ret)
3394 {
3395 int32 retval, Planet, itry;
3396 char ObjectName[AS_MAXCH], serr[AS_MAXCH], s[AS_MAXCH];
3397 double tjd0 = JDNDaysUTStart, tjd, dsynperiod, tjdmax, tadd;
3398 int32 MaxCountSynodicPeriod = MAX_COUNT_SYNPER;
3399 char *sevent[7] = {"", "morning first", "evening last", "evening first", "morning last", "acronychal rising", "acronychal setting"};
3400 if (dgeo[2] < SEI_ECL_GEOALT_MIN || dgeo[2] > SEI_ECL_GEOALT_MAX) {
3401 if (serr_ret != NULL)
3402 sprintf(serr_ret, "location for heliacal events must be between %.0f and %.0f m above sea\n", SEI_ECL_GEOALT_MIN, SEI_ECL_GEOALT_MAX);
3403 return ERR;
3404 }
3405 swi_set_tid_acc(JDNDaysUTStart, helflag, 0, serr);
3406 if (helflag & SE_HELFLAG_LONG_SEARCH)
3407 MaxCountSynodicPeriod = MAX_COUNT_SYNPER_MAX;
3408 /* if (helflag & SE_HELFLAG_SEARCH_1_PERIOD)
3409 MaxCountSynodicPeriod = 1; */
3410 if (serr_ret != NULL)
3411 *serr_ret = '\0';
3412 /* note, the fixed stars functions rewrite the star name. The input string
3413 may be too short, so we have to make sure we have enough space */
3414 strcpy_VBsafe(ObjectName, ObjectNameIn);
3415 tolower_string_star(ObjectName);
3416 default_heliacal_parameters(datm, dgeo, dobs, helflag);
3417 swe_set_topo(dgeo[0], dgeo[1], dgeo[2]);
3418 Planet = DeterObject(ObjectName);
3419 if (Planet == SE_SUN) {
3420 if (serr_ret != NULL) {
3421 strcpy(serr_ret, "the sun has no heliacal rising or setting\n");
3422 }
3423 return ERR;
3424 }
3425 /*
3426 * Moon events
3427 */
3428 if (Planet == SE_MOON) {
3429 if (TypeEvent == 1 || TypeEvent == 2) {
3430 if (serr_ret != NULL) {
3431 sprintf(serr_ret, "%s (event type %d) does not exist for the moon\n", sevent[TypeEvent], TypeEvent);
3432 }
3433 return ERR;
3434 }
3435 tjd = tjd0;
3436 retval = MoonEventJDut(tjd, dgeo, datm, dobs, TypeEvent, helflag, dret, serr);
3437 while (retval != -2 && *dret < tjd0) {
3438 tjd += 15;
3439 *serr = '\0';
3440 retval = MoonEventJDut(tjd, dgeo, datm, dobs, TypeEvent, helflag, dret, serr);
3441 }
3442 if (serr_ret != NULL && *serr != '\0')
3443 strcpy(serr_ret, serr);
3444 return retval;
3445 }
3446 /*
3447 * planets and fixed stars
3448 */
3449 if (!(helflag & SE_HELFLAG_AVKIND)) {
3450 if (Planet == -1 || Planet >= SE_MARS) {
3451 if (TypeEvent == 3 || TypeEvent == 4) {
3452 if (serr_ret != NULL) {
3453 if (Planet == -1)
3454 strcpy(s, ObjectName);
3455 else
3456 swe_get_planet_name(Planet, s);
3457 sprintf(serr_ret, "%s (event type %d) does not exist for %s\n", sevent[TypeEvent], TypeEvent, s);
3458 }
3459 return ERR;
3460 }
3461 }
3462 }
3463 /* arcus visionis method: set the TypeEvent for acronychal events */
3464 if (helflag & SE_HELFLAG_AVKIND) {
3465 if (Planet == -1 || Planet >= SE_MARS) {
3466 if (TypeEvent == SE_ACRONYCHAL_RISING)
3467 TypeEvent = 3;
3468 if (TypeEvent == SE_ACRONYCHAL_SETTING)
3469 TypeEvent = 4;
3470 }
3471 /* acronychal rising and setting (cosmic setting) are ill-defined.
3472 * We do not calculate them with the "visibility limit method" */
3473 } else if (1) {
3474 if (TypeEvent == SE_ACRONYCHAL_RISING || TypeEvent == SE_ACRONYCHAL_SETTING) {
3475 if (serr_ret != NULL) {
3476 if (Planet == -1)
3477 strcpy(s, ObjectName);
3478 else
3479 swe_get_planet_name(Planet, s);
3480 sprintf(serr_ret, "%s (event type %d) is not provided for %s\n", sevent[TypeEvent], TypeEvent, s);
3481 }
3482 return ERR;
3483 }
3484 }
3485 dsynperiod = get_synodic_period(Planet);
3486 tjdmax = tjd0 + dsynperiod * MaxCountSynodicPeriod;
3487 tadd = dsynperiod * 0.6;
3488 if (Planet == SE_MERCURY)
3489 tadd = 30;
3490 /*
3491 * this is the outer loop over n synodic periods
3492 */
3493 tjd = tjd0;
3494 retval = -2; /* indicates that another synodic period has to be done */
3495 for (itry = 0;
3496 tjd < tjdmax && retval == -2;
3497 itry++, tjd += tadd) {
3498 *serr = '\0';
3499 retval = heliacal_ut(tjd, dgeo, datm, dobs, ObjectName, TypeEvent, helflag, dret, serr);
3500 /* if resulting event date < start date for search (tjd0): retry starting
3501 * from half a period later. The event must be found now, unless there
3502 * is none, as is often the case with Mercury */
3503 while (retval != -2 && *dret < tjd0) {
3504 tjd += tadd;
3505 *serr = '\0';
3506 retval = heliacal_ut(tjd, dgeo, datm, dobs, ObjectName, TypeEvent, helflag, dret, serr);
3507 }
3508 }
3509 /*
3510 * no event was found within MaxCountSynodicPeriod, return error
3511 */
3512 if ((helflag & SE_HELFLAG_SEARCH_1_PERIOD) && (retval == -2 || dret[0] > tjd0 + dsynperiod * 1.5)) {
3513 strcpy(serr, "no heliacal date found within this synodic period");
3514 retval = -2;
3515 } else if (retval == -2) {
3516 sprintf(serr, "no heliacal date found within %d synodic periods", MaxCountSynodicPeriod);
3517 retval = ERR;
3518 }
3519 if (serr_ret != NULL && *serr != '\0')
3520 strcpy(serr_ret, serr);
3521 return retval;
3522 }
3523