1 /* NIGHTFALL Light Curve Synthesis Program */
2 /* Copyright (C) 1998 Rainer Wichmann */
3 /* */
4 /* This program is free software; you can redistribute it */
5 /* and/or modify */
6 /* it under the terms of the GNU General Public License as */
7 /* published by */
8 /* the Free Software Foundation; either version 2 of the License, or */
9 /* (at your option) any later version. */
10 /* */
11 /* This program is distributed in the hope that it will be useful, */
12 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
13 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
14 /* GNU General Public License for more details. */
15 /* */
16 /* You should have received a copy of the GNU General Public License */
17 /* along with this program; if not, write to the Free Software */
18 /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
19
20 #include <math.h>
21 #include <stdio.h>
22 #include <string.h>
23 #include <float.h>
24 #include "Light.h"
25
26
27
28 /****************************************************************************
29 @package nightfall
30 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
31 @version 1.0
32 @short Attach the throat of star a to star b
33 @tip This is to verify eclipse by 'own' star
34 @param (void)
35 @return (void)
36 @heading Light Curve
37 ****************************************************************************/
LightCopyThroat()38 void LightCopyThroat()
39 {
40 long j; /* loop counter */
41 long Nplus = 0; /* surface elements + throat */
42 long NP = 0; /* surface elements primary */
43 long NS = 0; /* surface elements secondary */
44
45 /* -------- attach Secondary throat to Primary ---------------------- */
46
47 Nplus = Binary[Secondary].N_PhiStep[StepsPerPi-1]
48 + Binary[Secondary].N_PhiStep[StepsPerPi-2];
49 NP = Binary[Primary].NumElem;
50 NS = Binary[Secondary].NumElem;
51
52
53 Binary[Primary].NumPlus = NP + Nplus;
54 if (Binary[Primary].NumPlus > (long) MaximumElements)
55 nf_error(_(errmsg[10]));
56
57 for (j = 0; j < Nplus; ++j) {
58
59 Surface[Primary][NP+j]
60 = Surface[Secondary][NS-Nplus+j];
61
62 /* --------------- INVERT X --------------------------------------- */
63
64 Surface[Primary][NP+j].lambda = (float)
65 (1.0 - Surface[Primary][NP+j].lambda);
66
67 Surface[Primary][NP+j].l =
68 (-Surface[Primary][NP+j].l);
69
70 }
71
72 /* ---------------- attach Primary throat to Secondary ---------------- */
73
74 Nplus = Binary[Primary].N_PhiStep[StepsPerPi-1]
75 + Binary[Primary].N_PhiStep[StepsPerPi-2];
76
77 Binary[Secondary].NumPlus = NS + Nplus;
78
79 if (Binary[Primary].NumPlus > (long) MaximumElements)
80 nf_error(_(errmsg[10]));
81
82 for (j = 0; j < Nplus; ++j) {
83
84 Surface[Secondary][NS+j]
85 = Surface[Primary][NP-Nplus+j];
86
87 /* --------------------- INVERT X --------------------------------- */
88
89 Surface[Secondary][NS+j].lambda = (float)
90 (1.0 - Surface[Secondary][NS+j].lambda);
91
92 Surface[Secondary][NS+j].l =
93 (-Surface[Secondary][NS+j].l);
94
95 }
96
97 return;
98 }
99
100 /****************************************************************************
101 @package nightfall
102 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
103 @version 1.0
104 @short Copy the throat back
105 @param (void)
106 @return (void)
107 @heading Light Curve
108 ****************************************************************************/
LightCopyThroatBack()109 void LightCopyThroatBack()
110 {
111 long j; /* loop counter */
112 long Nplus = 0; /* surface elements + throat */
113 long NP = 0; /* surface elements primary */
114 long NS = 0; /* surface elements secondary */
115
116 NP = Binary[Primary].NumElem;
117 NS = Binary[Secondary].NumElem;
118
119 /* ----------- copy back Primary throat to Primary ------------------- */
120
121 Nplus = Binary[Secondary].NumPlus - NS;
122
123 for (j = 0; j < Nplus; ++j) {
124
125 if(Surface[Secondary][NS+j].EclipseFlag == ON){
126
127 Surface[Primary][NP-Nplus+j].EclipseFlag
128 = Surface[Secondary][NS+j].EclipseFlag;
129
130 Surface[Primary][NP-Nplus+j].visibility
131 = Surface[Secondary][NS+j].visibility;
132
133 }
134
135 }
136
137 /* ------------- copy back Secondary throat to Secondary ------------- */
138
139 Nplus = Binary[Primary].NumPlus - NP;
140
141 for (j = 0; j < Nplus; ++j) {
142
143 if (Surface[Primary][NP+j].EclipseFlag == ON){
144
145 Surface[Secondary][NS-Nplus+j].EclipseFlag
146 = Surface[Primary][NP+j].EclipseFlag;
147
148 Surface[Secondary][NS-Nplus+j].visibility
149 = Surface[Primary][NP+j].visibility;
150
151 }
152 }
153
154 return;
155 }
156
157 /********************************************************************
158 @package nightfall
159 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
160 @version 1.0
161 @short Calculate various quantities for simple geometric tests
162 @param (void)
163 @return (int) Exit status
164 @heading Light Curve
165 ********************************************************************/
LightSetupTests()166 int LightSetupTests()
167 {
168
169 int testerr = 0; /* Exit status RootFind */
170
171 /* double Xl, h1, h2; */ /* Xl lagrange one */
172 double Xp1, Xp2; /* max. distance surface-center of star */
173 /* double PhiLagrange; */ /* tangent angle of roche lobes at L1 */
174 double CosPhiLagrange; /* tangent angle at Lagrange One (cos) */
175 double CosPhiOne = 0.0; /* perhaps more stringent if Roche lobe */
176 /* not filled */
177
178 /* ---------- simple geometric tests for eclipse possible ------ */
179
180
181 /* osculating cone angle at L1 */
182 /*
183 Xl = Binary[Primary].RLag1;
184 h1 = 1.0/(Xl*Xl*Xl);
185 h2 = 1.0/((1.0-Xl)*(1.0-Xl)*(1.0-Xl));
186 PhiLagrange = sqrt((2*h1 + 2*Mq*h2 + (1.0 - Mq))/(h1 + Mq*h2 - (1.0 - Mq)));
187 PhiLagrange = atan(PhiLagrange);
188 CosPhiLagrange = cos(PhiLagrange);
189 */
190
191 Xp1 = 1.0; Xp2 = 1.0;
192
193 /* calculate Xp == distance from center to surface towards L1 */
194
195 /* lambda, nu global; potential on x-axis */
196
197 lambda = 1.0; nu = 0.0;
198 RochePot = Binary[Primary].RochePot;
199 Mq = Binary[Primary].Mq;
200 if (Flags.asynchron1 == ON)
201 F = Binary[Primary].Fratio;
202 else
203 F = 1.0;
204
205
206 if (fabs(Binary[Primary].RocheFill - 1.0) <= FLT_EPSILON) {
207 Xp1 = Binary[Primary].RXCrit;
208 } else {
209 if (Flags.fill == OFF) {
210 Xp1 = RootFind(RocheSurface, 0.00001, Binary[Primary].RXCrit, 1.0e-8,
211 "LightSetupTests", &testerr);
212 if (testerr == 1) return(8);
213 } else {
214 Xp1 = Binary[Primary].LimRad;
215 }
216 }
217 Binary[Primary].Xp = Xp1;
218
219 RochePot = Binary[Secondary].RochePot;
220 Mq = Binary[Secondary].Mq;
221 if (Flags.asynchron2 == ON)
222 F = Binary[Secondary].Fratio;
223 else
224 F = 1.0;
225
226
227 if (fabs(Binary[Secondary].RocheFill - 1.0) <= FLT_EPSILON) {
228 Xp2 = Binary[Secondary].RXCrit;
229 } else {
230 if (Flags.fill == OFF) {
231 Xp2 = RootFind(RocheSurface, 0.00001,Binary[Secondary].RXCrit,
232 1.0e-8,
233 "LightSetupTests", &testerr );
234 if (testerr == 8) return(8);
235
236 } else {
237 Xp2 = Binary[Secondary].LimRad;
238 }
239 }
240 Binary[Secondary].Xp = Xp2;
241
242 if (Flags.debug[GEO] == ON)
243 {
244 printf("\nXp[P]: %6.2f Xp[S]: %6.2f\n",
245 Binary[Primary].Xp, Binary[Secondary].Xp);
246 }
247
248 /* ------ sanity check if F < 1.0 ------------------------------ */
249
250 if (
251 ( (Binary[Secondary].Fratio - 1.0) <= (-FLT_EPSILON)
252 && Flags.asynchron1 == ON )
253 || ( (Binary[Primary].Fratio - 1.0) <= (-FLT_EPSILON)
254 && Flags.asynchron2 == ON ) ) {
255 if ((Binary[Secondary].Xp + Binary[Primary].Xp) >= (1.0+DBL_EPSILON) ) {
256 WARNING(_("Fill Factor too large in Async Rotating Star --> Stars intersect"));
257 return(1);
258 }
259 }
260
261 if (Flags.fill == OFF) {
262 CosPhiOne = sqrt( MAX(0.0, (1.0 - SQR(Xp1 + Xp2))) );
263 } else {
264 CosPhiOne = 0.1;
265 }
266
267 /* osculating cone angle at L1 */
268
269 /* Chanan et al. (1972), ApJ 208, 512
270 */
271
272 if ( MIN(Binary[Primary].Mq, Binary[Secondary].Mq) > 0.4)
273 CosPhiLagrange = cos (DTOR*57.88);
274 else if ( MIN(Binary[Primary].Mq, Binary[Secondary].Mq) > 0.2)
275 CosPhiLagrange = cos (DTOR*59.00);
276 else if ( MIN(Binary[Primary].Mq, Binary[Secondary].Mq) > 0.1)
277 CosPhiLagrange = cos (DTOR*60.56);
278 else if ( MIN(Binary[Primary].Mq, Binary[Secondary].Mq) > 0.01)
279 CosPhiLagrange = cos (DTOR*67.17);
280 else
281 CosPhiLagrange = 0.0;
282
283 /* exact contact system
284 */
285 if ( (fabs(Binary[Primary].RocheFill - 1.0) <= FLT_EPSILON) &&
286 (fabs(Binary[Secondary].RocheFill - 1.0) <= FLT_EPSILON) ) {
287 CosPhiLagrange = CosPhiLagrange;
288 } else {
289 CosPhiLagrange = CosPhiOne;
290 }
291
292 Binary[Primary].CosPhiLagrange = CosPhiLagrange;
293 Binary[Secondary].CosPhiLagrange = CosPhiLagrange;
294 #ifdef HAVE_DISK
295 /*
296 * CosPhiLagrange for Disk; same formula as for CosPhiOne, but
297 * now with the outer disk radius.
298 */
299 Xp2 = Binary[Disk].Rout * Binary[Secondary].RLag1;
300 Binary[Disk].CosPhiLagrange = sqrt( MAX(0.0, (1.0 - SQR(Xp1 + Xp2))) );
301 #endif
302
303 return(0);
304 }
305
306 #ifdef HAVE_DISK
307
308 void SecondaryEclipsedByDisk (double CosPhase2, double SinPhase2,
309 double lzero, double mzero, double nzero);
310 void PrimaryEclipsedByDisk (double CosPhase, double SinPhase, double Phase,
311 double lzero, double mzero);
312 void DiskEclipsedByDisk (double CosPhase2, double SinPhase2,
313 double lzero, double mzero);
314
315
316 /********************************************************************
317 @package nightfall
318 @author Patrick Risse (risse@astro.uni-tuebingen.de)
319 @version 1.0
320 @short Eclipse Verification for a CB-System with a Disk
321 @param (int) j The phase index
322 @return (void)
323 @heading Light Curve
324 ********************************************************************/
LightCurveDisk(int j)325 void LightCurveDisk(int j)
326 {
327
328 long NElemD; /* # of surface elements */
329 long i; /* loop counter */
330 int eclipseS=0;
331 int eclipseP=0;
332 int eclipseD=0;
333 int count=0;
334 int invis=0;
335 int intersect = 0;
336 int testsec = 0;
337 int tested = 0;
338 int EclipseImpossible = 0;
339
340 int eclipsed = 0; /* test variable */
341
342 const int No = 0; /* test value */
343 const int Invisible = 10; /* test value */
344
345 int MinErr; /* exit status MinFinder */
346 double Phase, Phase2; /* orbital phase */
347
348 double lzero, mzero, nzero; /* LOS vector */
349 double lzero2, mzero2, nzero2; /* LOS vector */
350 double x_z, y_z, z_z;
351 double x_tl, y_tl, z_tl;
352 double Qi, Li, QLi; /* for sphere test */
353 double t_l, t_1, t_2, t_dist; /* intersection w/ sphere */
354
355 double tmin = 0; /* minimum along LOS */
356 double PotMin; /* minimum along LOS */
357 double SqrXp1; /* square of (center - lagrange one) */
358 double SqrPol1,SqrPol2; /* square of polar radius */
359 double Sini, Cosi; /* cos, sin of orbital inclination */
360 double SinPhase, CosPhase; /* cos, sin of phase (primary) */
361 double SinPhase2, CosPhase2; /* cos, sin of phase (secondary) */
362 double CosPhiLagrange; /* tangent angle at Lagrange One (cos) */
363
364 double Xp1, Xp2; /* max. distance surface-center of star */
365 double Pot_One, Pot_Two; /* surface potentials */
366
367 SurfaceElement *SurfPtrD; /* pointer to surface Disk */
368
369 double fratio1 = 1.0; /* async rotation */
370 double fratio2 = 1.0; /* async rotation */
371
372 double OmegaDot; /* velocity calculation */
373 double SecondaryMass; /* velocity calculation */
374 double SemiAxis; /* velocity calculation */
375 double KeplerPeriod; /* velocity calculation */
376 double G_SI = 6.6726e-11;/* gravitational constant */
377
378 /* COMMENTS IN ENGLISH PLEASE !!!!!!!!! MK */
379 /*double vec_e1[3];*/ /* Ebenenvektor 1 */
380 /*double vec_e2[3];*/ /* Ebenenvektor 2 */
381 /*double vec_p1[3];*/ /* Eckpunkt der Ebene */
382 /*double vec_en[3];*/ /* Normalenvektor der Ebene */
383 /*double vec_g1[3];*/ /* Geradenvektor */
384 /*double vec_p2[3];*/ /* Eckpunkt der Gerade */
385 /*double vec_dif_p2p1[3];*/ /* Differenze der beiten Ortsvektoren */
386 /*double det1,det2;*/ /* Determinanten */
387 /* double vec_intersection[3];*/ /* Schnittpunkt zwischen Gerade und Ebene */
388 /* double Distance[3]; */ /* Distanz zwischen Schnittpunkt und
389 Flaechenmittelpunkt */
390 /* double r; */ /* Faktor */
391 /* double rad_angle,angle; */ /* Angle between two vectors */
392
393 /* Recalculate the position of the Disk */
394 /* DivideDisk(Disk,j); */
395
396 if (Flags.asynchron1 == ON) fratio1 = Binary[Primary].Fratio;
397 if (Flags.asynchron2 == ON) fratio2 = Binary[Secondary].Fratio;
398
399
400 Sini = sin(Orbit.Inclination);
401 if (fabs(Sini) <= DBL_EPSILON)
402 Sini = 0.0;
403 Cosi = cos(Orbit.Inclination);
404 if (fabs(Cosi) <= DBL_EPSILON)
405 Cosi = 0.0;
406
407
408 Mq = Binary[Primary].Mq; /* Mq is global variable */
409
410 Pot_One = Binary[Primary].RochePot;
411 Pot_Two = Binary[Secondary].RochePot;
412 Xp1 = Binary[Primary].Xp;
413 Xp2 = Binary[Secondary].Xp;
414
415 CosPhiLagrange = Binary[Disk].CosPhiLagrange;
416
417 SqrXp1 = SQR(Xp1);
418
419 SqrPol1 = SQR(Binary[Primary].Radius); /* square of polar radius */
420 SqrPol2 = SQR(Binary[Secondary].Radius);
421
422 /* ------------------ eclipse test ------------------ */
423
424 /* to allow for generalization, we dont insist on */
425 /* same number of surface elements on both components */
426
427
428 NElemD = Binary[Disk].NumElem;
429
430 /* We start at 90 deg, Phase = -0.25, both stars fully visible */
431
432 if (Flags.elliptic == OFF) {
433
434 Phase = (PI/2 + (2*PI/PhaseSteps) * j );
435 FluxOut[j].Phase = (float) Phase;
436 Orbit.Phase = Phase;
437 Orbit.Nu = Phase;
438
439 } else {
440
441 Phase = Orbit.Phase;
442 /* Orbit.OmegaZero is mean anomaly at secondary eclipse */
443 FluxOut[j].Phase = (float) (Orbit.M + PI + Orbit.OmegaZero);
444 }
445
446 SinPhase = sin(Phase);
447 if (fabs(SinPhase) <= DBL_EPSILON)
448 SinPhase = 0.0;
449 CosPhase = cos(Phase);
450 if (fabs(CosPhase) <= DBL_EPSILON)
451 CosPhase = 0.0;
452
453 /* vector towards observer (LOS = Line Of Sight) */
454 lzero = Sini * CosPhase;
455 mzero = Sini * SinPhase;
456 nzero = Cosi;
457
458 /* Secondary turn by PI and mirror y->(-y) */
459 /* this is equivalent to mirror by x->(-x) */
460
461 Phase2 = (PI + Phase);
462 SinPhase2 = sin(Phase2);
463 if (fabs(SinPhase2) <= DBL_EPSILON)
464 SinPhase2 = 0.0;
465 CosPhase2 = cos(Phase2);
466 if (fabs(CosPhase2) <= DBL_EPSILON)
467 CosPhase2 = 0.0;
468
469 /* vector towards observer (LOS = Line Of Sight) */
470 lzero2 = Sini * CosPhase2;
471 mzero2 = Sini * SinPhase2;
472 nzero2 = Cosi;
473
474 SurfPtrD = Surface[Disk];
475
476
477 /* >>>>>>>>>>>>>>>>> the disk <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
478
479 /****** NEW CODE by MPR 29.04.2002 **********************/
480
481 /* Disk surrounds Secondary. Non-tilted display in Gnuplot ok. */
482
483 /* loop over surface elements */
484
485 for (i = 0; i < NElemD; ++i)
486 {
487 #if 0
488 if (fabs(SurfPtrD->lambda) <= DBL_EPSILON)
489 SurfPtrD->lambda = 0.0;
490 if (fabs(SurfPtrD->mu) <= DBL_EPSILON)
491 SurfPtrD->mu = 0.0;
492 if (fabs(SurfPtrD->nu) <= DBL_EPSILON)
493 SurfPtrD->nu = 0.0;
494 #endif
495
496 /* Need Kepler orbital period for disc element */
497
498 SecondaryMass = Orbit.TrueMass * Binary[Primary].Mq / ( 1.0 + Binary[Primary].Mq);
499 SemiAxis = Orbit.TrueDistance * SurfPtrD->rad;
500 KeplerPeriod = 2.0 * PI *
501 sqrt((SemiAxis * SemiAxis * SemiAxis) / (G_SI * SecondaryMass));
502
503 /* OmegaDot is rotation velocity in rad/s */
504
505 OmegaDot = 2.0 * PI / KeplerPeriod;
506
507 SurfPtrD->Velocity =
508 - OmegaDot * SemiAxis
509 * (SurfPtrD->mu * lzero2 + SurfPtrD->lambda * mzero2);
510
511 /* gamma = angle surface normal - LOS */
512
513 SurfPtrD->CosGamma = ( lzero2 * SurfPtrD->l
514 + mzero2 * (-SurfPtrD->m)
515 + nzero2 * SurfPtrD->n );
516
517 /*
518 if (i == 0) {
519 printf ("lzero2: %5.3f mzero2: %5.3f nzero2: %5.3f CosGamma: %5.3f\n",
520 lzero2, mzero2, nzero2,
521 SurfPtrD->CosGamma);
522 }
523 */
524
525 /*=========================================================== */
526
527 /* Default is visible */
528
529 eclipsed = No;
530
531 /* Start checking. */
532
533 /* 1. Does the surface normal point away from observer ?
534 * (Tested, works.)
535 */
536 if (SurfPtrD->CosGamma < DBL_EPSILON )
537 {
538 eclipsed = Invisible;
539 }
540
541 if (eclipsed == No)
542 {
543 /* 2. Is element eclipsed by secondary ?
544 * (Tested, works for non-inclined disk.)
545 *
546 * Test for an intersection with surface of Secondary.
547 * Both Secondary and disk are at (0,0,0). Phase is shifted
548 * by PI, and y->(-y).
549 */
550 x_z = SurfPtrD->lambda;
551 y_z = (-SurfPtrD->mu);
552 z_z = SurfPtrD->nu;
553 Qi = 2.0 * ( x_z*lzero2 + y_z*mzero2 + z_z*nzero2 );
554 Li = x_z*x_z + y_z*y_z + z_z*z_z - (Xp2)*(Xp2);
555 QLi = Qi*Qi - 4.0*Li;
556
557
558 /* If true there is an intersection with the
559 * sphere of the secondary star. Further tests
560 * are required
561 */
562 if (QLi <= 0.0)
563 {
564 /* no intersection */;
565 }
566 else
567 {
568
569 t_l = sqrt(QLi);
570 t_1 = (-Qi+t_l)/2;
571 t_2 = (-Qi-t_l)/2;
572
573 if (t_1 > 0.0 && t_2 > 0.0)
574 {
575 /* check whether distance is greater than polar radius */
576 /* -- maximum inscribed sphere */
577 t_l = (t_1 + t_2)/2.0;
578
579 /* t_l is Midpoint of ray intersecting Xp2-Sphere */
580 x_tl = x_z + t_l * lzero2;
581 y_tl = y_z + t_l * mzero2;
582 z_tl = z_z + t_l * nzero2;
583 t_dist = SQR(x_tl) + y_tl * y_tl + z_tl * z_tl;
584
585 if (t_dist <= (SqrPol2-DBL_EPSILON))
586 {
587 /* definitely eclipsed
588 */
589 eclipsed = EclipsedBySecondary;
590 }
591 else
592 {
593
594 /* Find Potential Minimum between t_1 and t_2
595 * final resort
596 * changed Binary[Pri].Mq to Binary[Sec].Mq Jul 6
597 */
598 MinErr = MinFinder(&Binary[Secondary].Mq,
599 &fratio2,
600 &t_1, &t_2,
601 &lzero, &mzero, &nzero,
602 &x_z, &y_z, &z_z,
603 &tmin, &PotMin);
604 if (MinErr == 1 && Flags.fill == OFF)
605 WARNING(_("LightCurve: No Minimum Found"));
606 PotMin = -PotMin;
607 SurfPtrD->LOS_Pot = (float) PotMin;
608
609 /* added condition tmin >= 0
610 */
611 if (PotMin >= Pot_Two && tmin >= 0)
612 {
613 eclipsed = EclipsedBySecondary;
614 }
615
616 }
617 } /* if (t_1 > 0.0 && t_2 > 0.0) */
618 } /* QLi > 0.0 */
619 }
620
621
622 if (eclipsed == No && CosPhase2 > (CosPhiLagrange-DBL_EPSILON))
623 {
624 ++testsec;
625
626 if (eclipsed == No)
627 {
628 /* 3. Is the element eclipsed by the primary?
629 * (Tested, works for non-inclined disk.)
630 */
631 x_z = SurfPtrD->lambda;
632 y_z = (-SurfPtrD->mu);
633 z_z = SurfPtrD->nu;
634
635 Qi = 2.0*( x_z*lzero2 + y_z*mzero2 + z_z*nzero2 - lzero2);
636 Li = 1.0 + x_z*x_z + y_z*y_z + z_z*z_z - 2*x_z - SqrXp1;
637 QLi = Qi*Qi - 4.0*Li;
638
639 if (QLi <= 0.0)
640 {
641 /* no intersection */;
642 }
643 else
644 {
645 t_l = sqrt(QLi); /* just to hold temporary result */
646 t_1 = (-Qi + t_l)/2.0;
647 t_2 = (-Qi - t_l)/2.0;
648
649 if (t_1 > 0.0 && t_2 > 0.0) /* intersection */
650 {
651 ++intersect;
652
653 /* check whether distance is greater than polar radius */
654 /* -- maximum inscribed sphere */
655 t_l = sqrt(QLi); /* just to hold temporary result */
656 t_1 = (-Qi + t_l)/2.0;
657 t_2 = (-Qi - t_l)/2.0;
658 t_l = (t_1 + t_2)/2.0;
659
660 /* t_l is Midpoint of ray intersecting Xp1-Sphere */
661 x_tl = x_z + t_l * lzero2;
662 y_tl = y_z + t_l * mzero2;
663 z_tl = z_z + t_l * nzero2;
664 t_dist = SQR(x_tl-1.0) + y_tl * y_tl + z_tl * z_tl;
665
666 /* we do not need the sqrt, compare the squares */
667 if (t_dist <= (SqrPol1-DBL_EPSILON))
668 {
669 /* definitely eclipsed
670 */
671 eclipsed = EclipsedByPrimary ;
672 }
673 else
674 {
675 ++tested;
676 /* Find Potential Minimum between t_1 and t_2 */
677 /* final resort */
678 /* changed Binary[Pri].Mq to Binary[Sec].Mq Jul 6 */
679 MinErr = MinFinder(&Binary[Primary].Mq,
680 &fratio1,
681 &t_1, &t_2,
682 &lzero2, &mzero2, &nzero2,
683 &x_z, &y_z, &z_z,
684 &tmin, &PotMin);
685 if (MinErr == 1)
686 WARNING(_("LightCurve: No Minimum Found"));
687 PotMin = -PotMin;
688 SurfPtrD->LOS_Pot = (float) PotMin;
689
690 /* added condition tmin >= 0 */
691 if (PotMin >= Pot_One && tmin >= 0)
692 {
693 eclipsed = EclipsedByPrimary;
694 }
695
696 } /* t_dist > (SqrPol1-DBL_EPSILON */
697 } /* if (t_1 > 0.0 && t_2 > 0.0) */
698 } /* QLi > 0.0) */
699 } /* (eclipsed < 1 )*/
700 }
701 else
702 {
703 ++EclipseImpossible;
704 }
705
706
707 if (eclipsed == EclipsedByPrimary)
708 eclipseP = eclipseP+1;
709 if (eclipsed == EclipsedByDisk)
710 eclipseD = eclipseD+1;
711 if (eclipsed == EclipsedBySecondary)
712 eclipseS = eclipseS+1;
713 if (eclipsed == Invisible)
714 invis = invis+1;
715 if (eclipsed == No)
716 {
717 count=count+1;
718 }
719
720
721 SurfPtrD->EclipseFlag = eclipsed;
722
723 if (SurfPtrD->EclipseFlag > 0 ) SurfPtrD->visibility = 0;
724 else SurfPtrD->visibility = 1;
725 if (eclipsed == 10) SurfPtrD->CosGamma = -1;
726 ++SurfPtrD;
727
728 } /* end loop over surface elements */
729
730 if (Flags.debug[GEO])
731 {
732 printf("=========================================\n");
733 printf("Step %02d\n", j);
734 printf("Total %d visible %d invisible %d eclipsed_tot %d\n",
735 count+invis+eclipseP+eclipseS+eclipseD,
736 count,invis,eclipseP+eclipseS+eclipseD);
737 printf("eclipsed: by Prim %d, by Sec %d, by Disk %d\n",
738 eclipseP,eclipseS,eclipseD);
739 printf("testsec : %d\n", testsec);
740 printf("intersect: %d\n", intersect);
741 printf("tested : %d\n", tested);
742 printf("imposs. : %d\n", EclipseImpossible);
743 }
744
745 count=0;
746 eclipseP=0;
747 eclipseS=0;
748 invis=0;
749 intersect=0;
750
751 SecondaryEclipsedByDisk (CosPhase2, SinPhase2, lzero2, mzero2, nzero2);
752 DiskEclipsedByDisk (CosPhase2, SinPhase2, lzero2, mzero2);
753
754 if (-CosPhase2 > (CosPhiLagrange-DBL_EPSILON))
755 PrimaryEclipsedByDisk (CosPhase, SinPhase, Phase, lzero, mzero);
756
757
758 return;
759
760 } /* end LightCurveDisk*/
761 #endif
762
763 /********************************************************************
764 @package nightfall
765 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
766 @version 1.0
767 @short Eclipse Verification
768 @param (int) j The phase index
769 @return (void)
770 @heading Light Curve
771 ********************************************************************/
LightCurve(int j)772 void LightCurve(int j)
773 {
774 const int primary_eclipsed = 1;
775 const int secondary_eclipsed = 2;
776 const int primary_total = 4;
777 const int secondary_total = 8;
778
779 long np_eclipsed = 0;
780 long ns_eclipsed = 0;
781 long np_visible = 0;
782 long ns_visible = 0;
783
784 static int contactflag = 0; /* changes at contacts */
785 int contactflag_old; /* to check for changes */
786
787 long NElem, NElemP, NElemS; /* # of surface elements */
788 long i; /* loop counter */
789 int eclipsed = 1; /* test variable */
790 const int Yes = 1; /* test value */
791 const int No = 0; /* test value */
792 const int Invisible = 10; /* test value */
793 int MinErr; /* exit status MinFinder */
794
795 double Phase, Phase2,M_Phase; /* orbital phase */
796 double OmegaDotPrimary; /* surface rotation async */
797 double OmegaDotSecondary; /* surface rotation async */
798 double lzero, mzero, nzero; /* LOS vector */
799 double lzero2, mzero2, nzero2; /* LOS vector */
800
801 double x_z, y_z, z_z;
802 double x_tl, y_tl, z_tl;
803 double Qi, Li, QLi; /* for sphere test */
804 double t_l, t_1, t_2, t_dist; /* intersection w/ sphere */
805
806 double tmin = 0; /* minimum along LOS */
807 double PotMin; /* minimum along LOS */
808 double SqrXp1, SqrXp2; /* square of (center - lagrange one) */
809 double SqrPol1,SqrPol2; /* square of polar radius */
810 double Sini, Cosi; /* cos, sin of orbital inclination */
811 double SinPhase, CosPhase; /* cos, sin of phase (primary) */
812 double SinPhase2, CosPhase2; /* cos, sin of phase (secondary) */
813 double Cosgamma; /* angle between LOS and surface normal */
814 double CosPhiLagrange; /* tangent angle at Lagrange One (cos) */
815
816 double Xp1, Xp2; /* max. distance surface-center of star */
817 double Pot_One, Pot_Two; /* surface potentials */
818
819 SurfaceElement *SurfPtrP; /* pointer to surface Primary */
820 SurfaceElement *SurfPtrS; /* pointer to surface Secondary */
821
822 double fratio1 = 1.0; /* async rotation */
823 double fratio2 = 1.0; /* async rotation */
824
825 if (Flags.asynchron1 == ON) fratio1 = Binary[Primary].Fratio;
826 if (Flags.asynchron2 == ON) fratio2 = Binary[Secondary].Fratio;
827
828 OmegaDotPrimary = fratio1 * 2.0 * PI / Orbit.TruePeriod;
829 OmegaDotSecondary = fratio2 * 2.0 * PI / Orbit.TruePeriod;
830
831 Sini = sin(Orbit.Inclination);
832 Cosi = cos(Orbit.Inclination);
833 NElem = Binary[Primary].NumElem;
834 Mq = Binary[Primary].Mq; /* Mq is global variable */
835
836 Pot_One = Binary[Primary].RochePot;
837 Pot_Two = Binary[Secondary].RochePot;
838 Xp1 = Binary[Primary].Xp;
839 Xp2 = Binary[Secondary].Xp;
840
841 /* CosPhiLagrange is same for Primary and Secondary */
842 CosPhiLagrange = Binary[Primary].CosPhiLagrange;
843
844 SqrXp1 = SQR(Xp1);
845 SqrXp2 = SQR(Xp2);
846 SqrPol1 = SQR(Binary[Primary].Radius); /* square of polar radius */
847 SqrPol2 = SQR(Binary[Secondary].Radius);
848
849 /* ------------------ eclipse test ------------------ */
850
851 /* to allow for generalization, we dont insist on */
852 /* same number of surface elements on both components */
853
854
855 NElem = IMAX(Binary[Primary].NumPlus,Binary[Secondary].NumPlus);
856 NElemP = IMAX(Binary[Primary].NumElem, Binary[Primary].NumPlus);
857 NElemS = IMAX(Binary[Secondary].NumElem, Binary[Secondary].NumPlus);
858
859 /*
860 NElem = IMAX(Binary[Primary].NumElem,Binary[Secondary].NumElem);
861 NElemP = Binary[Primary].NumElem;
862 NElemS = Binary[Secondary].NumElem;
863 */
864
865 /* We start at 90 deg, Phase = -0.25, both stars fully visible */
866
867 if (Flags.elliptic == OFF) {
868
869 Phase = (PI/2 + (2*PI/PhaseSteps) * j );
870 Orbit.Phase = Phase;
871 Orbit.Nu = Phase;
872 M_Phase = Phase;
873
874 } else {
875
876 Phase = Orbit.Phase;
877 /* Orbit.OmegaZero is mean anomaly at secondary eclipse */
878 M_Phase = Orbit.M + PI + Orbit.OmegaZero;
879 }
880
881 FluxOut[j].Phase = (float) M_Phase;
882
883 SinPhase = sin(Phase); CosPhase = cos(Phase);
884
885 /*
886 fprintf(stderr, _("#FIXME %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f (solar radius)\n"),
887 Orbit.Dist,
888 Orbit.Dist*Orbit.MinDist,
889 Orbit.TrueDistance / SOL_RADIUS,
890 Orbit.Dist * Binary[Primary].Radius * PERIDIST / SOL_RADIUS,
891 Binary[Primary].Radius * Orbit.TrueDistance / SOL_RADIUS,
892 Orbit.Dist * Binary[Primary].Xp * PERIDIST / SOL_RADIUS,
893 Binary[Primary].Xp * Orbit.TrueDistance / SOL_RADIUS,
894 Binary[Primary].RadMean * PERIDIST / SOL_RADIUS);
895 */
896
897 /* vector towards observer (LOS = Line Of Sight) */
898 lzero = Sini * CosPhase;
899 mzero = Sini * SinPhase;
900 nzero = Cosi;
901
902 /* Secondary turn by PI and mirror y->(-y) */
903 /* this is equivalent to mirror by x->(-x) */
904
905 Phase2 = (PI + Phase);
906 SinPhase2 = sin(Phase2); CosPhase2 = cos(Phase2);
907
908 /* vector towards observer (LOS = Line Of Sight) */
909 lzero2 = Sini * CosPhase2;
910 mzero2 = Sini * SinPhase2;
911 nzero2 = Cosi;
912
913 /* test whether eclipse of Primary definitely impossible */
914
915 /* this flag is (Flags.InEclipse1) required for fractional visibility
916 */
917 if ( CosPhase <= (CosPhiLagrange-DBL_EPSILON) ) {
918 Flags.InEclipse1 = OFF;
919 } else {
920 Flags.InEclipse1 = ON;
921 }
922
923
924 /* test whether eclipse of Secondary definitely impossible */
925
926 /* this flag (Flags.InEclipse2) is required for fractional visibility
927 */
928 if ( CosPhase2 <= (CosPhiLagrange-DBL_EPSILON) ) {
929 Flags.InEclipse2 = OFF;
930 } else {
931 Flags.InEclipse2 = ON;
932 }
933
934
935 SurfPtrP = Surface[Primary];
936 SurfPtrS = Surface[Secondary];
937
938 for (i = 0; i < NElem; ++i) { /* loop over surface elements */
939
940 /* START WITH PRIMARY */
941
942 if (i < NElemP) {
943
944 SurfPtrP->Velocity =
945 - OmegaDotPrimary * Orbit.TrueDistance * Orbit.Dist /* FIXME */
946 * (-SurfPtrP->mu * lzero
947 + SurfPtrP->lambda * mzero);
948
949
950 SurfPtrP->LOS_Pot = 0.0;
951 /* angle surface normal - LOS */
952 Cosgamma = lzero * SurfPtrP->l
953 + mzero * SurfPtrP->m
954 + nzero * SurfPtrP->n;
955
956 /* IF (Cosgamma > 0) THEN Visible Side, Test for Eclipse */
957 if (Cosgamma <= 0)
958 eclipsed = Invisible; /* invisible */
959 if (Cosgamma >= DBL_EPSILON){ /* Visible Side */
960
961 eclipsed = Yes;
962 ++np_visible;
963
964 /* test whether eclipse definitely impossible */
965 if (CosPhase <= (CosPhiLagrange-DBL_EPSILON)) eclipsed = No;
966
967 if (Flags.fill == ON) {
968 if (i > Binary[Primary].NumElem) eclipsed = Yes;
969 }
970
971 /* IF (eclipsed == Yes) MORE TESTING */
972 if (eclipsed == Yes) {
973
974 /* intersect with circle of rad Xp2 centered on second */
975 /* -- minimum sphere containing the star */
976 x_z = SurfPtrP->lambda;
977 y_z = SurfPtrP->mu;
978 z_z = SurfPtrP->nu;
979 Qi = 2.0*( x_z*lzero + y_z*mzero + z_z*nzero - lzero);
980 Li = 1.0 + x_z*x_z + y_z*y_z + z_z*z_z - 2*x_z - SqrXp2;
981 QLi = Qi*Qi - 4.0*Li;
982 if (QLi <= 0.0) { /* No intersection */
983 eclipsed = No;
984 } else { /* Intersection, next test */
985
986 /* check whether distance is greater than polar radius */
987 /* -- maximum inscribed sphere */
988 t_l = sqrt(QLi); /* just to hold temporary result */
989 t_1 = (-Qi + t_l)/2.0;
990 t_2 = (-Qi - t_l)/2.0;
991 t_l = (t_1 + t_2)/2.0;
992 /* t_l is Midpoint of ray intersecting Xp2-Sphere */
993 x_tl = x_z + t_l * lzero;
994 y_tl = y_z + t_l * mzero;
995 z_tl = z_z + t_l * nzero;
996 t_dist = SQR(x_tl-1.0) + y_tl * y_tl + z_tl * z_tl;
997
998 if (Flags.fill == ON) {
999 if (i > Binary[Primary].NumElem) {
1000 /* rw 06.02.2007 */
1001 /* eclipsed = Yes; */ /* superfluous */
1002 /* incorrect, will test against wrong star
1003 * t_1 = MAX(t_1, t_2);
1004 * t_2 = t_1 + 1.0;
1005 */
1006
1007 /* next test would yield incorrect result, thus jump
1008 * to the MinFinder() routine
1009 */
1010 t_dist = SqrPol2+1.;
1011 }
1012 }
1013
1014 if (t_dist <= (SqrPol2-DBL_EPSILON)) {
1015
1016 eclipsed = Yes; /* definitely eclipsed */
1017
1018 } else {
1019
1020
1021 /* Find Potential Minimum between t_1 and t_2 */
1022 /* final resort */
1023 /* changed Binary[Pri].Mq to Binary[Sec].Mq Jul 6 */
1024 MinErr = MinFinder(&Binary[Secondary].Mq,
1025 &fratio2,
1026 &t_1, &t_2,
1027 &lzero, &mzero, &nzero,
1028 &x_z, &y_z, &z_z,
1029 &tmin, &PotMin);
1030 if (MinErr == 1 && Flags.fill == OFF)
1031 WARNING(_("LightCurve: No Minimum Found"));
1032 PotMin = -PotMin;
1033 SurfPtrP->LOS_Pot = (float) PotMin;
1034
1035 /* added condition tmin >= 0 */
1036 if (PotMin >= Pot_Two && tmin >= 0) {
1037 eclipsed = Yes;
1038 } else {
1039 eclipsed = No;
1040 }
1041
1042 }
1043 }
1044 }
1045
1046 if (eclipsed == Yes) ++np_eclipsed;
1047 SurfPtrP->EclipseFlag = eclipsed;
1048 SurfPtrP->CosGamma = (float) Cosgamma;
1049 if (eclipsed > 0) SurfPtrP->visibility = 0;
1050 else SurfPtrP->visibility = 1;
1051
1052 } else {
1053 /* invisible side */
1054 SurfPtrP->EclipseFlag = eclipsed;
1055 SurfPtrP->CosGamma = -1;
1056 if (eclipsed > 0) SurfPtrP->visibility = 0;
1057 else SurfPtrP->visibility = 1;
1058 }
1059 ++SurfPtrP;
1060 }
1061 /* END PRIMARY */
1062
1063 /* HERE SECONDARY */
1064 if (i < NElemS) {
1065
1066 SurfPtrS->Velocity =
1067 - OmegaDotSecondary * Orbit.TrueDistance * Orbit.Dist /* FIXME */
1068 * (SurfPtrS->mu * lzero2
1069 + SurfPtrS->lambda * mzero2);
1070
1071 SurfPtrS->LOS_Pot = 0.0;
1072 /* angle surface normal - LOS */
1073 /* dont forget to mirror in y->(-y) */
1074 Cosgamma = lzero2 * SurfPtrS->l
1075 + mzero2 * (-SurfPtrS->m)
1076 + nzero2 * SurfPtrS->n;
1077
1078 /* IF (Cosgamma > 0) THEN Visible Side, Test for Eclipse */
1079 if (Cosgamma <= 0)
1080 eclipsed = Invisible; /* invisible */
1081 if (Cosgamma >= DBL_EPSILON){ /* Visible Side */
1082
1083 eclipsed = Yes;
1084 ++ns_visible;
1085
1086 /* test whether eclipse definitely impossible */
1087 if (CosPhase2 <= (CosPhiLagrange-DBL_EPSILON)) eclipsed = No;
1088
1089 if (Flags.fill == ON) {
1090 if (i > Binary[Secondary].NumElem) eclipsed = Yes;
1091 }
1092
1093 /* IF (eclipsed == Yes) MORE TESTING */
1094 if (eclipsed == Yes) {
1095
1096 /* intersect with circle of rad Xp1 centered on prim */
1097 /* -- minimum sphere containing the star */
1098 x_z = SurfPtrS->lambda;
1099 y_z = (-SurfPtrS->mu);
1100 z_z = SurfPtrS->nu;
1101 Qi = 2.0*( x_z*lzero2 + y_z*mzero2 + z_z*nzero2 - lzero2);
1102 Li = 1.0 + x_z*x_z + y_z * y_z + z_z*z_z - 2*x_z - SqrXp1;
1103 QLi = Qi*Qi - 4.0*Li;
1104 if (QLi <= 0.0) { /* No intersection */
1105 eclipsed = No;
1106 } else { /* Intersection, next test */
1107
1108 /* check whether distance is greater than polar radius */
1109 /* -- maximum inscribed sphere */
1110 t_l = sqrt(QLi); /* just to hold temporary result */
1111 t_1 = (-Qi + t_l)/2.0;
1112 t_2 = (-Qi - t_l)/2.0;
1113 t_l = (t_1 + t_2)/2.0;
1114 /* t_l is Midpoint of ray intersecting Xp2-Sphere */
1115 x_tl = x_z + t_l * lzero2;
1116 y_tl = y_z + t_l * mzero2;
1117 z_tl = z_z + t_l * nzero2;
1118 t_dist = SQR(x_tl-1.0) + y_tl*y_tl + z_tl*z_tl;
1119
1120 if (Flags.fill == ON) {
1121 if (i > Binary[Secondary].NumElem) {
1122 /* rw 06.02.2007 */
1123 /* eclipsed = Yes; */ /* superfluous */
1124 /* incorrect, will test against wrong star
1125 * t_1 = MAX(t_1, t_2);
1126 * t_2 = t_1 + 1.0;
1127 */
1128
1129 /* next test would yield incorrect result, thus jump
1130 * to the MinFinder() routine
1131 */
1132 t_dist = SqrPol1+1.;
1133 }
1134 }
1135
1136 /* we do not need the sqrt, compare the squares */
1137 if (t_dist <= (SqrPol1-DBL_EPSILON)) {
1138
1139 eclipsed = Yes; /* definitely eclipsed */
1140
1141 } else {
1142
1143 /* Find Potential Minimum between t_1 and t_2 */
1144 /* final resort */
1145 /* changed Binary[Sec].Mq to Binary[Pri].Mq Jul 6 */
1146 MinErr = MinFinder(&Binary[Primary].Mq,
1147 &fratio1,
1148 &t_1, &t_2,
1149 &lzero2, &mzero2, &nzero2,
1150 &x_z, &y_z, &z_z,
1151 &tmin, &PotMin);
1152 if (MinErr == 1) WARNING(_("LightCurve: No Minimum Found"));
1153 PotMin = -PotMin;
1154 SurfPtrS->LOS_Pot = (float) PotMin;
1155
1156 /* added condition tmin >= 0 */
1157 if (PotMin >= (Pot_One-DBL_EPSILON) && tmin >= 0) {
1158 eclipsed = Yes;
1159 } else {
1160 eclipsed = No;
1161 }
1162
1163
1164 }
1165 }
1166 }
1167
1168 if (eclipsed == Yes) ++ns_eclipsed;
1169 SurfPtrS->EclipseFlag = eclipsed;
1170 SurfPtrS->CosGamma = (float) Cosgamma;
1171 if (eclipsed > 0) SurfPtrS->visibility = 0;
1172 else SurfPtrS->visibility = 1;
1173
1174 } else {
1175 SurfPtrS->EclipseFlag = eclipsed;
1176 SurfPtrS->CosGamma = -1;
1177 if (eclipsed > 0) SurfPtrS->visibility = 0;
1178 else SurfPtrS->visibility = 1;
1179 }
1180 ++SurfPtrS;
1181 }
1182 /* END SECONDARY */
1183
1184
1185 } /* end loop over surface elements */
1186
1187 M_Phase = Orbit.M + PI + Orbit.OmegaZero;
1188
1189 contactflag_old = contactflag;
1190
1191 if (ns_eclipsed > 0)
1192 contactflag |= secondary_eclipsed;
1193 else
1194 contactflag &= ~secondary_eclipsed;
1195
1196 if (ns_eclipsed == ns_visible)
1197 contactflag |= secondary_total;
1198 else
1199 contactflag &= ~secondary_total;
1200
1201 if (np_eclipsed > 0)
1202 contactflag |= primary_eclipsed;
1203 else
1204 contactflag &= ~primary_eclipsed;
1205
1206 if (np_eclipsed == np_visible)
1207 contactflag |= primary_total;
1208 else
1209 contactflag &= ~primary_total;
1210
1211 if ((contactflag & secondary_eclipsed) != (contactflag_old & secondary_eclipsed)) {
1212 if (contactflag & secondary_eclipsed) {
1213 if (Flags.debug[VERBOSE] == ON)
1214 printf(_(" Secondary eclipse first contact at phase: %8.6f\n"), Phase);
1215 Orbit.Contact[1] = M_Phase;
1216 }
1217 else {
1218 if (Flags.debug[VERBOSE] == ON)
1219 printf(_(" Secondary eclipse fourth contact at phase: %8.6f\n"), Phase);
1220 Orbit.Contact[2] = M_Phase;
1221 }
1222 }
1223
1224 if ((contactflag & secondary_total) != (contactflag_old & secondary_total)) {
1225 if (contactflag & secondary_total) {
1226 if (Flags.debug[VERBOSE] == ON)
1227 printf(_(" Secondary eclipse second contact at phase: %8.6f\n"), Phase);
1228 Orbit.Contact[5] = M_Phase;
1229 }
1230 else {
1231 if (Flags.debug[VERBOSE] == ON)
1232 printf(_(" Secondary eclipse third contact at phase: %8.6f\n"), Phase);
1233 Orbit.Contact[6] = M_Phase;
1234 }
1235 }
1236
1237 if ((contactflag & primary_eclipsed) != (contactflag_old & primary_eclipsed)) {
1238 if (contactflag & primary_eclipsed) {
1239 if (Flags.debug[VERBOSE] == ON)
1240 printf(_(" Primary eclipse first contact at phase: %8.6f\n"), Phase);
1241 Orbit.Contact[3] = M_Phase;
1242 }
1243 else {
1244 if (Flags.debug[VERBOSE] == ON)
1245 printf(_(" Primary eclipse fourth contact at phase: %8.6f\n"), Phase);
1246 Orbit.Contact[4] = M_Phase;
1247 }
1248 }
1249
1250 if ((contactflag & primary_total) != (contactflag_old & primary_total)) {
1251 if (contactflag & primary_total) {
1252 if (Flags.debug[VERBOSE] == ON)
1253 printf(_(" Primary eclipse second contact at phase: %8.6f\n"), Phase);
1254 Orbit.Contact[7] = M_Phase;
1255 }
1256 else {
1257 if (Flags.debug[VERBOSE] == ON)
1258 printf(_(" Primary eclipse third contact at phase: %8.6f\n"), Phase);
1259 Orbit.Contact[8] = M_Phase;
1260 }
1261 }
1262
1263 return;
1264 } /* end LightCurve */
1265
1266 /* From an article by Mark Lowne / mark marklowne com
1267 * This allows a branch-free "if (x >= 0) {...}".
1268 * No speedup though on the "if (CosGamma >= 0) ..."
1269 */
ispositive(float x)1270 inline float ispositive(float x)
1271 {
1272 union {
1273 float f;
1274 int i;
1275 } v;
1276 v.f = x;
1277 v.i = ((~v.i) >> 8) & 0x3f800000;
1278
1279 return v.f;
1280 }
1281
1282 /********************************************************************
1283 @package nightfall
1284 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1285 @version 1.0
1286 @short Flux computation
1287 @param (int) Comp The stellar component
1288 @param (int) Phase The phase index
1289 @return (void)
1290 @heading Light Curve
1291 ********************************************************************/
LightFlux(int Comp,int Phase)1292 void LightFlux(int Comp, int Phase)
1293 {
1294 register long i; /* loop counter */
1295 double Cosgamma; /* LOS angle */
1296 double Common; /* common part */
1297 double Weight; /* weight factor */
1298 double SumOfVelocity; /* normalization */
1299 double SumOfWeight; /* normalization */
1300 double OneMinCosgamma; /* 1.0 - Cosgamma */
1301 SurfaceElement *SurfPtr; /* pointer to surface */
1302 long nelem; /* loop limit */
1303 PhotoOut *OFlux; /* pointer to outflux */
1304
1305 /* ------------------ flux computing --------------- */
1306
1307 /* We start at 90 deg, Phase = 0.25, both stars fully visible */
1308
1309 /* loops over passbands are unrolled, all if's in loop eliminated */
1310
1311 SumOfVelocity = 0.0;
1312 SumOfWeight = 0.0;
1313 nelem = Binary[Comp].NumElem;
1314 OFlux = &FluxOut[Phase];
1315 SurfPtr = Surface[Comp];
1316
1317 if (Flags.limb == 0) {
1318
1319 for (i = 0; i < nelem; ++i) {
1320
1321 /* loop over surface elements */
1322 /* add up Flux of visible elements */
1323 /* IF NOT eclipsed THEN */
1324
1325 if (SurfPtr->CosGamma >= 0.0) {
1326 Cosgamma = SurfPtr->CosGamma;
1327 Common = SurfPtr->visibility
1328 * Cosgamma;
1329
1330 OneMinCosgamma = (1.0 - Cosgamma);
1331
1332 Weight = SurfPtr->f_[0] * Common
1333 * (1.0 - Limb[Comp][0][0]*OneMinCosgamma);
1334 OFlux->Mag[0] = OFlux->Mag[0] + Weight;
1335
1336 Weight = SurfPtr->f_[1] * Common
1337 * (1.0 - Limb[Comp][1][0]*OneMinCosgamma);
1338 OFlux->Mag[1] = OFlux->Mag[1] + Weight;
1339
1340 Weight = SurfPtr->f_[2] * Common
1341 * (1.0 - Limb[Comp][2][0]*OneMinCosgamma);
1342 OFlux->Mag[2] = OFlux->Mag[2] + Weight;
1343
1344 SumOfVelocity = SumOfVelocity
1345 + Weight * SurfPtr->Velocity;
1346 SumOfWeight = SumOfWeight + Weight;
1347
1348 Weight = SurfPtr->f_[3] * Common
1349 * (1.0 - Limb[Comp][3][0]*OneMinCosgamma);
1350 OFlux->Mag[3] = OFlux->Mag[3] + Weight;
1351
1352 Weight = SurfPtr->f_[4] * Common
1353 * (1.0 - Limb[Comp][4][0]*OneMinCosgamma);
1354 OFlux->Mag[4] = OFlux->Mag[4] + Weight;
1355
1356 Weight = SurfPtr->f_[5] * Common
1357 * (1.0 - Limb[Comp][5][0]*OneMinCosgamma);
1358 OFlux->Mag[5] = OFlux->Mag[5] + Weight;
1359
1360 Weight = SurfPtr->f_[6] * Common
1361 * (1.0 - Limb[Comp][6][0]*OneMinCosgamma);
1362 OFlux->Mag[6] = OFlux->Mag[6] + Weight;
1363
1364 Weight = SurfPtr->f_[7] * Common
1365 * (1.0 - Limb[Comp][7][0]*OneMinCosgamma);
1366 OFlux->Mag[7] = OFlux->Mag[7] + Weight;
1367
1368 Weight = SurfPtr->f_[8] * Common
1369 * (1.0 - Limb[Comp][8][0]*OneMinCosgamma);
1370 OFlux->Mag[8] = OFlux->Mag[8] + Weight;
1371
1372 Weight = SurfPtr->f_[9] * Common
1373 * (1.0 - Limb[Comp][9][0]*OneMinCosgamma);
1374 OFlux->Mag[9] = OFlux->Mag[9] + Weight;
1375
1376 Weight = SurfPtr->f_[10] * Common
1377 * (1.0 - Limb[Comp][10][0]*OneMinCosgamma);
1378 OFlux->Mag[10] = OFlux->Mag[10] + Weight;
1379
1380 Weight = SurfPtr->f_[11] * Common
1381 * (1.0 - Limb[Comp][10][0]*OneMinCosgamma);
1382 OFlux->Mag[11] = OFlux->Mag[11] + Weight;
1383
1384 }
1385 ++SurfPtr;
1386 }
1387 }
1388
1389 else if ( Flags.limb == 1) {
1390
1391 double OneMinCosgamma2;
1392
1393 for (i = 0; i < nelem; ++i) {
1394
1395 /* loop over surface elements */
1396 /* add up Flux of visible elements */
1397 /* IF NOT eclipsed THEN */
1398
1399 if (SurfPtr->CosGamma >= 0.0) {
1400
1401 Cosgamma = SurfPtr->CosGamma;
1402 Common = SurfPtr->visibility
1403 * Cosgamma;
1404 OneMinCosgamma = (1.0 - Cosgamma);
1405 OneMinCosgamma2 = OneMinCosgamma*OneMinCosgamma;
1406
1407 Weight = SurfPtr->f_[0] * Common
1408 *(1.0 - Limb[Comp][0][0]*OneMinCosgamma
1409 - Limb[Comp][0][1]*OneMinCosgamma2);
1410 OFlux->Mag[0] = OFlux->Mag[0] + Weight;
1411 Weight = SurfPtr->f_[1] * Common
1412 *(1.0 - Limb[Comp][1][0]*OneMinCosgamma
1413 - Limb[Comp][1][1]*OneMinCosgamma2);
1414 OFlux->Mag[1] = OFlux->Mag[1] + Weight;
1415 Weight = SurfPtr->f_[2] * Common
1416 *(1.0 - Limb[Comp][2][0]*OneMinCosgamma
1417 - Limb[Comp][2][1]*OneMinCosgamma2);
1418 OFlux->Mag[2] = OFlux->Mag[2] + Weight;
1419
1420 SumOfVelocity = SumOfVelocity
1421 + Weight * SurfPtr->Velocity;
1422 SumOfWeight = SumOfWeight + Weight;
1423
1424 Weight = SurfPtr->f_[3] * Common
1425 *(1.0 - Limb[Comp][3][0]*OneMinCosgamma
1426 - Limb[Comp][3][1]*OneMinCosgamma2);
1427 OFlux->Mag[3] = OFlux->Mag[3] + Weight;
1428 Weight = SurfPtr->f_[4] * Common
1429 *(1.0 - Limb[Comp][4][0]*OneMinCosgamma
1430 - Limb[Comp][4][1]*OneMinCosgamma2);
1431 OFlux->Mag[4] = OFlux->Mag[4] + Weight;
1432 Weight = SurfPtr->f_[5] * Common
1433 *(1.0 - Limb[Comp][5][0]*OneMinCosgamma
1434 - Limb[Comp][5][1]*OneMinCosgamma2);
1435 OFlux->Mag[5] = OFlux->Mag[5] + Weight;
1436 Weight = SurfPtr->f_[6] * Common
1437 *(1.0 - Limb[Comp][6][0]*OneMinCosgamma
1438 - Limb[Comp][6][1]*OneMinCosgamma2);
1439 OFlux->Mag[6] = OFlux->Mag[6] + Weight;
1440 Weight = SurfPtr->f_[7] * Common
1441 *(1.0 - Limb[Comp][7][0]*OneMinCosgamma
1442 - Limb[Comp][7][1]*OneMinCosgamma2);
1443 OFlux->Mag[7] = OFlux->Mag[7] + Weight;
1444 Weight = SurfPtr->f_[8] * Common
1445 *(1.0 - Limb[Comp][8][0]*OneMinCosgamma
1446 - Limb[Comp][8][1]*OneMinCosgamma2);
1447 OFlux->Mag[8] = OFlux->Mag[8] + Weight;
1448 Weight = SurfPtr->f_[9] * Common
1449 *(1.0 - Limb[Comp][9][0]*OneMinCosgamma
1450 - Limb[Comp][9][1]*OneMinCosgamma2);
1451 OFlux->Mag[9] = OFlux->Mag[9] + Weight;
1452 Weight = SurfPtr->f_[10] * Common
1453 *(1.0 - Limb[Comp][10][0]*OneMinCosgamma
1454 - Limb[Comp][10][1]*OneMinCosgamma2);
1455 OFlux->Mag[10] = OFlux->Mag[10] + Weight;
1456 Weight = SurfPtr->f_[11] * Common
1457 *(1.0 - Limb[Comp][11][0]*OneMinCosgamma
1458 - Limb[Comp][11][1]*OneMinCosgamma2);
1459 OFlux->Mag[11] = OFlux->Mag[11] + Weight;
1460
1461 }
1462 ++SurfPtr;
1463 }
1464 }
1465 else if ( Flags.limb == 2) {
1466
1467 double sqrtCosgamma;
1468
1469 for (i = 0; i < nelem; ++i) {
1470
1471 /* loop over surface elements */
1472 /* add up Flux of visible elements */
1473 /* IF NOT eclipsed THEN */
1474
1475 if (SurfPtr->CosGamma >= 0.0) {
1476
1477 Cosgamma = SurfPtr->CosGamma;
1478 Common = SurfPtr->visibility
1479 * Cosgamma;
1480 OneMinCosgamma = (1.0 - Cosgamma);
1481 sqrtCosgamma = 1.0 - sqrt(Cosgamma);
1482
1483 Weight = SurfPtr->f_[0] * Common
1484 *(1.0 - Limb[Comp][0][0]*OneMinCosgamma
1485 - Limb[Comp][0][1]*sqrtCosgamma);
1486 OFlux->Mag[0] = OFlux->Mag[0] + Weight;
1487 Weight = SurfPtr->f_[1] * Common
1488 *(1.0 - Limb[Comp][1][0]*OneMinCosgamma
1489 - Limb[Comp][1][1]*sqrtCosgamma);
1490 OFlux->Mag[1] = OFlux->Mag[1] + Weight;
1491 Weight = SurfPtr->f_[2] * Common
1492 *(1.0 - Limb[Comp][2][0]*OneMinCosgamma
1493 - Limb[Comp][2][1]*sqrtCosgamma);
1494 OFlux->Mag[2] = OFlux->Mag[2] + Weight;
1495
1496 SumOfVelocity = SumOfVelocity
1497 + Weight * SurfPtr->Velocity;
1498 SumOfWeight = SumOfWeight + Weight;
1499
1500 Weight = SurfPtr->f_[3] * Common
1501 *(1.0 - Limb[Comp][3][0]*OneMinCosgamma
1502 - Limb[Comp][3][1]*sqrtCosgamma);
1503 OFlux->Mag[3] = OFlux->Mag[3] + Weight;
1504 Weight = SurfPtr->f_[4] * Common
1505 *(1.0 - Limb[Comp][4][0]*OneMinCosgamma
1506 - Limb[Comp][4][1]*sqrtCosgamma);
1507 OFlux->Mag[4] = OFlux->Mag[4] + Weight;
1508 Weight = SurfPtr->f_[5] * Common
1509 *(1.0 - Limb[Comp][5][0]*OneMinCosgamma
1510 - Limb[Comp][5][1]*sqrtCosgamma);
1511 OFlux->Mag[5] = OFlux->Mag[5] + Weight;
1512 Weight = SurfPtr->f_[6] * Common
1513 *(1.0 - Limb[Comp][6][0]*OneMinCosgamma
1514 - Limb[Comp][6][1]*sqrtCosgamma);
1515 OFlux->Mag[6] = OFlux->Mag[6] + Weight;
1516 Weight = SurfPtr->f_[7] * Common
1517 *(1.0 - Limb[Comp][7][0]*OneMinCosgamma
1518 - Limb[Comp][7][1]*sqrtCosgamma);
1519 OFlux->Mag[7] = OFlux->Mag[7] + Weight;
1520 Weight = SurfPtr->f_[8] * Common
1521 *(1.0 - Limb[Comp][8][0]*OneMinCosgamma
1522 - Limb[Comp][8][1]*sqrtCosgamma);
1523 OFlux->Mag[8] = OFlux->Mag[8] + Weight;
1524 Weight = SurfPtr->f_[9] * Common
1525 *(1.0 - Limb[Comp][9][0]*OneMinCosgamma
1526 - Limb[Comp][9][1]*sqrtCosgamma);
1527 OFlux->Mag[9] = OFlux->Mag[9] + Weight;
1528 Weight = SurfPtr->f_[10] * Common
1529 *(1.0 - Limb[Comp][10][0]*OneMinCosgamma
1530 - Limb[Comp][10][1]*sqrtCosgamma);
1531 OFlux->Mag[10] = OFlux->Mag[10] + Weight;
1532 Weight = SurfPtr->f_[11] * Common
1533 *(1.0 - Limb[Comp][11][0]*OneMinCosgamma
1534 - Limb[Comp][11][1]*sqrtCosgamma);
1535 OFlux->Mag[11] = OFlux->Mag[11] + Weight;
1536
1537 }
1538 ++SurfPtr;
1539 }
1540 }
1541 else if ( Flags.limb == 3) {
1542
1543 double sqrtCosgamma;
1544 double Cosgamma2;
1545 double Cosgamma32;
1546
1547 for (i = 0; i < nelem; ++i) {
1548
1549 /* loop over surface elements */
1550 /* add up Flux of visible elements */
1551 /* IF NOT eclipsed THEN */
1552
1553 if (SurfPtr->CosGamma >= 0.0) {
1554
1555 Cosgamma = SurfPtr->CosGamma;
1556 Common = SurfPtr->visibility
1557 * Cosgamma;
1558
1559 sqrtCosgamma = 1.0 - sqrt(Cosgamma);
1560 OneMinCosgamma = 1.0 - Cosgamma;
1561 Cosgamma32 = 1.0 - Cosgamma*sqrt(Cosgamma);
1562 Cosgamma2 = 1.0 - Cosgamma*Cosgamma;
1563
1564 Weight = SurfPtr->f_[0] * Common
1565 *(1.0 - Limb[Comp][0][0]*sqrtCosgamma - Limb[Comp][0][1]*OneMinCosgamma
1566 - Limb[Comp][0][2]*Cosgamma32 - Limb[Comp][0][3]*Cosgamma2);
1567 OFlux->Mag[0] = OFlux->Mag[0] + Weight;
1568 Weight = SurfPtr->f_[1] * Common
1569 *(1.0 - Limb[Comp][1][0]*sqrtCosgamma - Limb[Comp][1][1]*OneMinCosgamma
1570 - Limb[Comp][1][2]*Cosgamma32 - Limb[Comp][1][3]*Cosgamma2);
1571 OFlux->Mag[1] = OFlux->Mag[1] + Weight;
1572 Weight = SurfPtr->f_[2] * Common
1573 *(1.0 - Limb[Comp][2][0]*sqrtCosgamma - Limb[Comp][2][1]*OneMinCosgamma
1574 - Limb[Comp][2][2]*Cosgamma32 - Limb[Comp][2][3]*Cosgamma2);
1575 OFlux->Mag[2] = OFlux->Mag[2] + Weight;
1576
1577 SumOfVelocity = SumOfVelocity
1578 + Weight * SurfPtr->Velocity;
1579 SumOfWeight = SumOfWeight + Weight;
1580
1581 Weight = SurfPtr->f_[3] * Common
1582 *(1.0 - Limb[Comp][3][0]*sqrtCosgamma - Limb[Comp][3][1]*OneMinCosgamma
1583 - Limb[Comp][3][2]*Cosgamma32 - Limb[Comp][3][3]*Cosgamma2);
1584 OFlux->Mag[3] = OFlux->Mag[3] + Weight;
1585 Weight = SurfPtr->f_[4] * Common
1586 *(1.0 - Limb[Comp][4][0]*sqrtCosgamma - Limb[Comp][4][1]*OneMinCosgamma
1587 - Limb[Comp][4][2]*Cosgamma32 - Limb[Comp][4][3]*Cosgamma2);
1588 OFlux->Mag[4] = OFlux->Mag[4] + Weight;
1589 Weight = SurfPtr->f_[5] * Common
1590 *(1.0 - Limb[Comp][5][0]*sqrtCosgamma - Limb[Comp][5][1]*OneMinCosgamma
1591 - Limb[Comp][5][2]*Cosgamma32 - Limb[Comp][5][3]*Cosgamma2);
1592 OFlux->Mag[5] = OFlux->Mag[5] + Weight;
1593 Weight = SurfPtr->f_[6] * Common
1594 *(1.0 - Limb[Comp][6][0]*sqrtCosgamma - Limb[Comp][6][1]*OneMinCosgamma
1595 - Limb[Comp][6][2]*Cosgamma32 - Limb[Comp][6][3]*Cosgamma2);
1596 OFlux->Mag[6] = OFlux->Mag[6] + Weight;
1597 Weight = SurfPtr->f_[7] * Common
1598 *(1.0 - Limb[Comp][7][0]*sqrtCosgamma - Limb[Comp][7][1]*OneMinCosgamma
1599 - Limb[Comp][7][2]*Cosgamma32 - Limb[Comp][7][3]*Cosgamma2);
1600 OFlux->Mag[7] = OFlux->Mag[7] + Weight;
1601 Weight = SurfPtr->f_[8] * Common
1602 *(1.0 - Limb[Comp][8][0]*sqrtCosgamma - Limb[Comp][8][1]*OneMinCosgamma
1603 - Limb[Comp][8][2]*Cosgamma32 - Limb[Comp][8][3]*Cosgamma2);
1604 OFlux->Mag[8] = OFlux->Mag[8] + Weight;
1605 Weight = SurfPtr->f_[9] * Common
1606 *(1.0 - Limb[Comp][9][0]*sqrtCosgamma - Limb[Comp][9][1]*OneMinCosgamma
1607 - Limb[Comp][9][2]*Cosgamma32 - Limb[Comp][9][3]*Cosgamma2);
1608 OFlux->Mag[9] = OFlux->Mag[9] + Weight;
1609 Weight = SurfPtr->f_[10] * Common
1610 *(1.0 - Limb[Comp][10][0]*sqrtCosgamma - Limb[Comp][10][1]*OneMinCosgamma
1611 - Limb[Comp][10][2]*Cosgamma32 - Limb[Comp][10][3]*Cosgamma2);
1612 OFlux->Mag[10] = OFlux->Mag[10] + Weight;
1613 Weight = SurfPtr->f_[11] * Common
1614 *(1.0 - Limb[Comp][11][0]*sqrtCosgamma - Limb[Comp][11][1]*OneMinCosgamma
1615 - Limb[Comp][11][2]*Cosgamma32 - Limb[Comp][11][3]*Cosgamma2);
1616 OFlux->Mag[11] = OFlux->Mag[11] + Weight;
1617
1618 }
1619 ++SurfPtr;
1620 }
1621 }
1622
1623 /* ------------ normalize velocity ---------------------- */
1624
1625 if (SumOfWeight >= DBL_EPSILON) {
1626 FluxOut[Phase].D_RV[Comp] = (float) (SumOfVelocity / SumOfWeight);
1627 } else {
1628 FluxOut[Phase].D_RV[Comp] = 0.0;
1629 }
1630
1631 if (Flags.elliptic == ON)
1632 FluxOut[Phase].D_RV[Comp] = (-1.0) * FluxOut[Phase].D_RV[Comp];
1633
1634 FluxOut[Phase].RV[Comp] =
1635 FluxOut[Phase].D_RV[Comp] + FluxOut[Phase].RV[Comp];
1636
1637 return;
1638 }
1639
1640 /****************************************************************************
1641 @package nightfall
1642 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1643 @version 1.0
1644 @short Find minimum of the potential along some line-of-sight vector.
1645 @long Uses Brents algorithm, which is a combination of golden section
1646 search and parabolic interpolation.
1647 Adopted from procedure 'localmin' in: Richard Brent, Algorithms for
1648 minimization without derivatives, Prentice-Hall, Inc. (1973)
1649 @tip Input is a parametric form of LOS vector + brackets.
1650 @param (double*) q Mass ratio
1651 @param (double*) ff Async rotation
1652 @param (double*) t1 Minimum bracket
1653 @param (double*) t3 Minimum bracket
1654 @param (double*) l0 x-component LOS vector
1655 @param (double*) m0 y-component LOS vector
1656 @param (double*) n0 z-component LOS vector
1657 @param (double*) x0 x-component of point on LOS vector
1658 @param (double*) y0 y-component of point on LOS vector
1659 @param (double*) z0 z-component of point on LOS vector
1660 @param (double*) tmin The position of the minimum
1661 @param (double*) Cmin The LOS potential minimum
1662 @return (int) The error status
1663 @heading Light Curve
1664 ****************************************************************************/
MinFinder(double * q,double * ff,double * t1,double * t3,double * l0,double * m0,double * n0,double * x0,double * y0,double * z0,double * tmin,double * Cmin)1665 int MinFinder(double *q, double *ff, double *t1, double *t3,
1666 double *l0, double *m0, double *n0,
1667 double *x0, double *y0, double *z0,
1668 double *tmin, /*@out@*/ double *Cmin)
1669 {
1670 register int i; /* loop variable */
1671 double t2; /* middle value */
1672 double tmp; /* temp var */
1673 double step, stepstore; /* stepsizes */
1674 double Fe = 0.0, Ff = 0.0; /* function values */
1675 double Fg = 0.0, Fh = 0.0; /* function values */
1676 double e = 0.0, f, g, h; /* positions (parametric) */
1677 double delta1, delta2; /* test convergence */
1678 double Start, End, Middle; /* brackets */
1679 double nextstep = 0.0; /* stepsize */
1680 double x, y, z; /* position in cartesian coordinates */
1681 double P1 = 0.0, P2, P3; /* fit parabola */
1682 double precise = 0.01; /* convergence criterion */
1683 double x2, y2, z2; /* squared values */
1684
1685 if ( (*t1) <= ((*t3)-DBL_EPSILON) ) {
1686 Start = (*t1);
1687 End = (*t3);
1688 } else {
1689 Start = (*t3);
1690 End = (*t1);
1691 }
1692
1693 step = 0.0;
1694
1695 /* we have a parametric description of a vector X(h) */
1696 /* we want to find min(func) = min(X(h)) as function of h */
1697
1698 t2 = 0.62*(*t1) + 0.38*(*t3);
1699
1700 /* initial value of X(h) */
1701 x = 1.0 - ((*x0) + (*l0)*(t2));
1702 y = -(*y0) - (*m0)*(t2);
1703 z = (*z0) + (*n0)*(t2);
1704
1705 x2 = x*x;
1706 y2 = y*y;
1707 z2 = z*z;
1708
1709 /* NEW jul 6,98 -- change of coordinate system, calculate */
1710 /* in system of eclipsing star */
1711 /* (i.e. x = 1.0 - x; y = -y; z = z;) */
1712
1713 /* t2 is the initial guess for the minimum between t1 and t3 */
1714 h = g = f = t2;
1715
1716 /* initial value of F */
1717 tmp = 1.0-x;
1718 tmp = sqrt(tmp*tmp + y2 + z2 );
1719 Ff = Fg = Fh = -( 1.0/sqrt(x2 + y2 + z2)
1720 + (*q) * (1.0/tmp - x)
1721 + ( x2 + y2 ) * (*ff)*(*ff) * (0.5*( (*q) + 1.0)));
1722
1723 /* we do a maximum of MAXITER iterations - convergence usually */
1724 /* after less than ten iterations */
1725 /*@i@*/ for (i = MAXITER; i; --i) {
1726
1727 Middle = 0.5 * (Start + End);
1728
1729 delta1 = precise * fabs(h) + DBL_EPSILON;
1730 delta2 = 2.0 * delta1;
1731
1732 /* test for convergence */
1733 /* RETURN if converged */
1734
1735 if (fabs(h-Middle) <= (delta2 - 0.5 * (End-Start) ) ) {
1736 *tmin = h;
1737 *Cmin = Fh;
1738 return(0);
1739 }
1740
1741 if (fabs(step) >= delta1) {
1742
1743 /* fit parabola */
1744 P3 = (h-g)*(Fh-Ff);
1745 P2 = (h-f)*(Fh-Fg);
1746 P1 = (h-f)*P2 - (h-g)*P3;
1747 P2 = 2.0 * (P2 - P3);
1748
1749 if (P2 >= DBL_EPSILON) {
1750 P1 = -P1;
1751 } else {
1752 P2 = -P2;
1753 }
1754
1755 stepstore = step;
1756 step = nextstep;
1757
1758 if ( (fabs(P1) >= fabs(0.5 * P2 * stepstore))
1759 || (P1 <= P2*(Start-h))
1760 || (P1 >= P2*(End-h))) {
1761
1762 /* a golden-section step */
1763 if (h >= Middle) {
1764 step = (Start-h);
1765 } else {
1766 step = (End-h);
1767 }
1768 nextstep = 0.381966 * step;
1769
1770 } else {
1771
1772 /* do a parabolic interpolation step */
1773 nextstep = P1/P2;
1774 e = h + nextstep;
1775
1776 /* f must not be evaluated too close to Start or End */
1777 if ((e-Start) <= delta2 || (End-e) <= delta2) {
1778 if ((Middle-h) <= (-DBL_EPSILON)) {
1779 nextstep = nextstep-delta1;
1780 } else {
1781 nextstep = nextstep+delta1;
1782 }
1783 }
1784
1785 }
1786
1787 } else {
1788
1789 /* a golden-section step */
1790 if (h >= Middle) {
1791 step = (Start-h);
1792 } else {
1793 step = (End-h);
1794 }
1795 nextstep = 0.381966 * step;
1796
1797 }
1798
1799 if (fabs(nextstep) >= delta1) {
1800 e = h + nextstep;
1801 } else {
1802 /* f must not be evaluated too close */
1803 if ( nextstep <= (-DBL_EPSILON) ) {
1804 e = h - delta1;
1805 } else {
1806 e = h + delta1;
1807 }
1808 }
1809
1810 /* e ist the current best guess for the minimum */
1811 /* calculate X(e) */
1812 x = 1.0 - ((*x0) + (*l0) * e);
1813 y = - (*y0) - (*m0) * e;
1814 z = (*z0) + (*n0) * e;
1815
1816 x2 = x*x;
1817 y2 = y*y;
1818 z2 = z*z;
1819
1820 /* NEW jul 6 -- change of coordinate system, calculate */
1821 /* in system of eclipsing star */
1822 /* i.e x = 1.0 - x; y = -y; z = z; */
1823
1824 /* calculate f(X(e)) */
1825 tmp = 1.0-x;
1826 tmp = sqrt(tmp*tmp + y2 + z2 );
1827 Fe = -( 1.0/sqrt(x2 + y2 + z2)
1828 + (*q) * (1.0/tmp - x)
1829 + ( x2 + y2 )* (*ff)*(*ff) * (0.5*( (*q) + 1.0)));
1830
1831 /* update */
1832 if (Fe <= Fh) {
1833 if ( e >= h ) Start = h;
1834 else End = h;
1835 f = g; g = h; h = e;
1836 Ff = Fg; Fg = Fh; Fh = Fe;
1837 } else {
1838 if ( e <= (h-DBL_EPSILON) ) Start = e;
1839 else End = e;
1840 if ( (Fe <= Fg) || (fabs (g - x) <= FLT_EPSILON) ) {
1841 f = g; g = e;
1842 Ff = Fg; Fg = Fe;
1843 }
1844 else if ( (Fe <= Ff) || (fabs (f- g) <= FLT_EPSILON)
1845 || (fabs (f - h) <= FLT_EPSILON) ) {
1846 f = e; Ff = Fe;
1847 }
1848 }
1849
1850 /* end of main loop */
1851 }
1852
1853 /* only reach here if no convergence */
1854 *tmin = h;
1855 *Cmin = Fh;
1856 return(1);
1857 }
1858
1859 #ifdef HAVE_DISK
1860 /********************************************************************
1861 @package nightfall
1862 @author Rainer Wichmann
1863 @version 1.69
1864 @short Eclipse checking
1865 @param (double) x coordinate
1866 @param (double) y coordinate
1867 @param (double) R radius
1868 @return (double) angle on circle
1869 @heading Light Curve
1870 ********************************************************************/
toAngle(double x,double y,double R)1871 static double toAngle (double x, double y, double R)
1872 {
1873 double angle;
1874
1875 x = x / R;
1876
1877 if (y < 0)
1878 angle = PI - acos(x);
1879 else
1880 angle = PI + acos(x);
1881
1882 return angle;
1883 }
1884
1885 /********************************************************************
1886 @package nightfall
1887 @author Rainer Wichmann
1888 @version 1.70
1889 @short Provide disk heights for eclipse testing
1890 @param (double*) Height of inner edge of the disk
1891 @param (double*) Height of outer edge of the disk
1892 @param (double*) Height of the puffed hot spot
1893 @return (void)
1894 @heading Light Curve
1895 ********************************************************************/
diskHeight(double * Height_in,double * Height_out,double * Height_puff)1896 void diskHeight (double * Height_in, double * Height_out, double * Height_puff)
1897 {
1898 double diskExp = 1.0; /* Exponent for r dependency of
1899 scale height */
1900 double diskProp = 1.0; /* Proportionality factor for disk
1901 height */
1902 double diskAdd = 1.0; /* Additive constant for disk
1903 height */
1904
1905 double Rout = Binary[Disk].RLag1 * Binary[Disk].Rout;
1906 double Rin = Binary[Disk].RLag1 * Binary[Disk].Rin;
1907
1908
1909 if (Flags.tdisk == 0)
1910 {
1911 diskAdd = Binary[Disk].Thick;
1912 diskProp = Binary[Disk].HR;
1913 diskExp = 1.0;
1914 }
1915 else
1916 {
1917 if (Flags.tdisk == 1) /* isothermal */
1918 diskExp = 1.5;
1919 else
1920 diskExp = (9.0/8.0); /* reprocessing */
1921
1922 diskAdd = 0.0;
1923 diskProp = (Binary[Disk].Thick / pow(Rin, diskExp));
1924 }
1925
1926 /* Height of the Disk on the outer side
1927 */
1928 *Height_out = (diskAdd + diskProp * pow(Rout, diskExp)) / 2.0;
1929
1930 /* Height of the Disk on the inner side
1931 */
1932 *Height_in = (diskAdd + diskProp * pow(Rin, diskExp)) / 2.0;
1933
1934 if (Binary[Disk].puffHS > FLT_EPSILON)
1935 *Height_puff = *Height_out + Binary[Disk].puffHS;
1936 else
1937 *Height_puff = *Height_out;
1938
1939 return;
1940 }
1941
1942 /* This macro yields correct results only if a > 0 and b in [0, 2*PI).
1943 */
1944 #define DIST_CIRCLE(a, b) fabs(PI - fabs(PI - fmod((a) - (b), (2.0 * PI) ) ))
1945
obscuredByPuff(double norm_dist,double height,double height_puff)1946 int obscuredByPuff(double norm_dist, double height, double height_puff)
1947 {
1948 if (norm_dist <= 1 && height <= height_puff && height >= (-height_puff))
1949 {
1950 double comp = puffHeight(Binary[Disk].puffHS, norm_dist);
1951 if (height <= comp && height >= (-comp))
1952 return 1;
1953 }
1954 return 0;
1955 }
1956
1957 /********************************************************************
1958 @package nightfall
1959 @author Rainer Wichmann
1960 @version 1.55
1961 @short Test eclipse of secondary by surrounding disk
1962 @param (double) CosPhase2 Cosine of phase
1963 @param (double) SinPhase2 Sine of phase
1964 @return (void)
1965 @heading Light Curve
1966 ********************************************************************/
SecondaryEclipsedByDisk(double CosPhase2,double SinPhase2,double lzero,double mzero,double nzero)1967 void SecondaryEclipsedByDisk (double CosPhase2, double SinPhase2,
1968 double lzero, double mzero, double nzero)
1969 {
1970 /* Radii of outer and inner cylinder between which
1971 * the disk is situated.
1972 */
1973 double Rout = Binary[Disk].RLag1 * Binary[Disk].Rout;
1974 double Rin = Binary[Disk].RLag1 * Binary[Disk].Rin;
1975 double R0;
1976
1977 /* Height of the Disk on the outer side
1978 */
1979 double HeightIn;
1980 double HeightOut;
1981 double HeightPuff;
1982
1983 double x_z;
1984 double y_z;
1985 double z_z;
1986
1987 double A;
1988 double B;
1989 double A2B;
1990 double temp;
1991
1992 SurfaceElement *SurfPtrS = Surface[Secondary];
1993
1994 unsigned long NElemS = Binary[Secondary].NumElem;
1995 unsigned long countS = 0;
1996
1997 double Tani = tan(Orbit.Inclination);
1998 double Sini = sin(Orbit.Inclination);;
1999 double Zin, Zout;
2000 double t_1;
2001 double t_2;
2002 double t_1_z;
2003 double t_2_z;
2004
2005 double hs_long = dtor_2pi(Binary[Disk].longitudeHS); /* longitude of hot spot */
2006 double hs_ext = ((Binary[Disk].extentHS / 2.0) * DTOR); /* extent in phi */
2007
2008 /* View from above
2009 */
2010 if (Tani < FLT_EPSILON)
2011 return;
2012
2013 diskHeight (&HeightIn, &HeightOut, &HeightPuff);
2014
2015 while (countS < NElemS)
2016 {
2017 /* Invisible anyway
2018 */
2019 if (SurfPtrS->visibility == 0)
2020 goto next;
2021
2022 /* Above disk (and above puff) => visible
2023 */
2024 if (SurfPtrS->nu > HeightPuff)
2025 goto next;
2026
2027 x_z = SurfPtrS->lambda;
2028 y_z = (-SurfPtrS->mu);
2029 z_z = SurfPtrS->nu;
2030
2031 /* Compute intersection with inner cylinder
2032 */
2033 A = x_z * CosPhase2 + y_z * SinPhase2;
2034 R0 = (x_z * x_z) + (y_z * y_z);
2035 B = R0 - (Rin * Rin);
2036 A2B = (A*A) - B;
2037
2038 if (A2B <= 0)
2039 goto next;
2040
2041 /* z-coordinate of intersection
2042 */
2043 temp = sqrt(A2B);
2044
2045 t_1 = ((-A) - temp) / Tani;
2046 t_2 = ((-A) + temp) / Tani;
2047
2048 Zin = MAX(t_1,t_2) + z_z;
2049
2050 /* Compute intersection with outer cylinder
2051 */
2052 B = R0 - (Rout * Rout);
2053 A2B = (A*A) - B;
2054
2055 if (A2B <= 0)
2056 goto next;
2057
2058 /* Z-coordinate of intersection. Here, t_1, t_2 are
2059 * t_1 * sin(i), t_2 * sin(i) for the solution
2060 * of the equation in x-y plane.
2061 */
2062 temp = sqrt(A2B);
2063
2064 t_1 = ((-A) - temp);
2065 t_2 = ((-A) + temp);
2066
2067 t_1_z = t_1 / Tani;
2068 t_2_z = t_2 / Tani;
2069
2070 Zout = MAX(t_1_z,t_2_z) + z_z;
2071
2072 if ( (Zout < (-HeightPuff)) || (Zin > HeightIn && Zout > HeightPuff) )
2073 goto next;
2074
2075 if ((Zin > HeightIn && Zout > HeightOut && Zout <= HeightPuff) ||
2076 (Zout < (-HeightOut) && Zout >= (-HeightPuff)))
2077 {
2078 if (Sini > FLT_EPSILON)
2079 {
2080 double a, hs_dist, t_xy;
2081
2082 if (fabs(t_1) < fabs(t_2))
2083 {
2084 t_xy = t_1 / Sini;
2085 }
2086 else
2087 {
2088 t_xy = t_2 / Sini;
2089 }
2090
2091 a = toAngle((x_z + t_xy*lzero), (y_z + t_xy*mzero), Rout);
2092
2093 hs_dist = DIST_CIRCLE(a, hs_long);
2094
2095 if (hs_dist <= hs_ext)
2096 {
2097 if (obscuredByPuff(fabs(hs_dist/hs_ext), Zout, HeightPuff))
2098 {
2099 goto set_eclipsed;
2100 }
2101 }
2102 }
2103 goto next;
2104 }
2105
2106 set_eclipsed:
2107 SurfPtrS->EclipseFlag = 1;
2108 SurfPtrS->visibility = 0;
2109
2110 next:
2111 ++SurfPtrS;
2112 ++countS;
2113 }
2114
2115 return;
2116 }
2117
2118 /********************************************************************
2119 @package nightfall
2120 @author Rainer Wichmann
2121 @version 1.55
2122 @short Test eclipse of disk by itself
2123 @param (double) CosPhase2 Cosine of phase
2124 @param (double) SinPhase2 Sine of phase
2125 @return (void)
2126 @heading Light Curve
2127 ********************************************************************/
DiskEclipsedByDisk(double CosPhase2,double SinPhase2,double lzero,double mzero)2128 void DiskEclipsedByDisk (double CosPhase2, double SinPhase2, double lzero, double mzero)
2129 {
2130 /* Radii of outer and inner cylinder between which
2131 * the disk is situated.
2132 */
2133 double Rout = Binary[Disk].RLag1 * Binary[Disk].Rout;
2134 double Rin = Binary[Disk].RLag1 * Binary[Disk].Rin;
2135 double R0;
2136
2137 /* Height of the Disk on the outer side
2138 */
2139 double HeightIn;
2140 double HeightOut;
2141 double HeightPuff;
2142
2143 double x_z;
2144 double y_z;
2145 double z_z;
2146
2147 double A;
2148 double B;
2149 double A2B;
2150
2151 SurfaceElement *SurfPtr = Surface[Disk];
2152
2153 unsigned long NElem = Binary[Disk].NumElem;
2154 unsigned long count = 0;
2155
2156 double Tani = tan(Orbit.Inclination);
2157 double Sini = sin(Orbit.Inclination);
2158 double Zin, Zout;
2159 double t_1;
2160 double t_2;
2161 double t_1_z;
2162 double t_2_z;
2163 double temp;
2164
2165 double hs_long = dtor_2pi(Binary[Disk].longitudeHS); /* longitude of hot spot */
2166 double hs_ext = ((Binary[Disk].extentHS / 2.0) * DTOR); /* extent in phi */
2167
2168 /* View from above
2169 */
2170 if (Tani < FLT_EPSILON)
2171 return;
2172
2173 diskHeight (&HeightIn, &HeightOut, &HeightPuff);
2174
2175 while (count < NElem)
2176 {
2177 /* Invisible anyway
2178 */
2179 if (SurfPtr->visibility == 0)
2180 goto next;
2181
2182 /* Above disk => visible
2183 */
2184 if (SurfPtr->nu > HeightPuff)
2185 goto next;
2186
2187 /* Outer side => never eclipsed by disk
2188 */
2189 if (SurfPtr->AreaType == OUT_RECTANGLE)
2190 goto next;
2191
2192 x_z = SurfPtr->lambda;
2193 y_z = (-SurfPtr->mu);
2194 z_z = SurfPtr->nu;
2195
2196 A = x_z * CosPhase2 + y_z * SinPhase2;
2197 R0 = (x_z * x_z) + (y_z * y_z);
2198
2199 if (SurfPtr->AreaType == IN_RECTANGLE)
2200 {
2201 /* Compute intersection with inner cylinder
2202 */
2203 B = R0 - (Rin * Rin);
2204 A2B = (A*A) - B;
2205
2206 if (A2B <= 0)
2207 goto outer;
2208
2209 /* z-coordinate of intersection
2210 */
2211 temp = sqrt(A2B);
2212
2213 t_1 = ((-A) - temp) / Tani;
2214 t_2 = ((-A) + temp) / Tani;
2215
2216 Zin = MAX(t_1,t_2) + z_z;
2217
2218 if (Zin < HeightIn && Zin > (-HeightIn))
2219 goto set_eclipsed;
2220 }
2221
2222 outer:
2223
2224 /* Compute intersection with outer cylinder
2225 */
2226 B = R0 - (Rout * Rout);
2227 A2B = (A*A) - B;
2228
2229 if (A2B <= 0)
2230 goto next;
2231
2232 /* z-coordinate of intersection
2233 */
2234 temp = sqrt(A2B);
2235
2236 t_1 = ((-A) - temp);
2237 t_2 = ((-A) + temp);
2238
2239 t_1_z = t_1 / Tani;
2240 t_2_z = t_2 / Tani;
2241
2242 Zout = MAX(t_1_z,t_2_z) + z_z;
2243
2244 if (Zout <= HeightOut)
2245 goto set_eclipsed;
2246
2247 if (Zout > HeightPuff)
2248 goto next;
2249
2250 /* Between HeightOut and HeightPuff, locate
2251 * intersecting point.
2252 */
2253 if (Sini > FLT_EPSILON)
2254 {
2255 double a, hs_dist, t_xy;
2256
2257 if (t_1_z > t_2_z)
2258 {
2259 t_xy = t_1 / Sini;
2260 }
2261 else
2262 {
2263 t_xy = t_2 / Sini;
2264 }
2265
2266 a = toAngle((x_z + t_xy*lzero), (y_z + t_xy*mzero), Rout);
2267
2268 hs_dist = DIST_CIRCLE(a, hs_long);
2269
2270 if (hs_dist <= hs_ext)
2271 {
2272 if (obscuredByPuff(fabs(hs_dist/hs_ext), Zout, HeightPuff))
2273 {
2274 goto set_eclipsed;
2275 }
2276 }
2277 }
2278
2279 goto next;
2280
2281
2282 set_eclipsed:
2283 SurfPtr->EclipseFlag = 1;
2284 SurfPtr->visibility = 0;
2285
2286 next:
2287 ++SurfPtr;
2288 ++count;
2289 }
2290
2291 return;
2292 }
2293
2294
2295
2296 /********************************************************************
2297 @package nightfall
2298 @author Rainer Wichmann
2299 @version 1.55
2300 @short Test eclipse of primary by disk surrounding secondary
2301 @param (double) CosPhase2 Cosine of phase
2302 @param (double) SinPhase2 Sine of phase
2303 @return (void)
2304 @heading Light Curve
2305 ********************************************************************/
PrimaryEclipsedByDisk(double CosPhase,double SinPhase,double Phase,double lzero,double mzero)2306 void PrimaryEclipsedByDisk (double CosPhase, double SinPhase, double Phase,
2307 double lzero, double mzero)
2308 {
2309 /* Radii of outer and inner cylinder between which
2310 * the disk is situated.
2311 */
2312 double Rout = Binary[Disk].RLag1 * Binary[Disk].Rout;
2313 double Rin = Binary[Disk].RLag1 * Binary[Disk].Rin;
2314
2315 /* Height of the Disk on the outer side
2316 */
2317 double HeightIn;
2318 double HeightOut;
2319 double HeightPuff;
2320
2321 double x_z;
2322 double y_z;
2323 double z_z;
2324
2325 double A;
2326 double B;
2327 double A2B;
2328 double R0;
2329
2330 SurfaceElement *SurfPtr = Surface[Primary];
2331
2332 unsigned long NElem = Binary[Primary].NumElem;
2333 unsigned long count = 0;
2334
2335 double Tani = tan(Orbit.Inclination);
2336 double Sini = sin(Orbit.Inclination);
2337 double Zmin;
2338 double Zmax;
2339 double t_1;
2340 double t_2;
2341 double t_1_z;
2342 double t_2_z;
2343 double temp;
2344
2345 double hs_long = dtor_2pi(Binary[Disk].longitudeHS); /* longitude of hot spot */
2346 double hs_ext = ((Binary[Disk].extentHS / 2.0) * DTOR); /* extent in phi */
2347
2348 /* View from above
2349 */
2350 if (Tani < FLT_EPSILON)
2351 return;
2352
2353 diskHeight (&HeightIn, &HeightOut, &HeightPuff);
2354
2355 while (count < NElem)
2356 {
2357 /* Invisible anyway
2358 */
2359 if (SurfPtr->visibility == 0)
2360 goto next;
2361
2362 /* Above disk => visible
2363 */
2364 if (SurfPtr->nu > HeightPuff)
2365 goto next;
2366
2367 x_z = SurfPtr->lambda;
2368 y_z = SurfPtr->mu;
2369 z_z = SurfPtr->nu;
2370
2371 /* Compute intersection with outer cylinder
2372 */
2373 R0 = 1.0 + (x_z * x_z) + (y_z * y_z) - (2.0 * x_z);
2374
2375 A = x_z * CosPhase + y_z * SinPhase - CosPhase;
2376 B = R0 - (Rout * Rout);
2377 A2B = (A*A) - B;
2378
2379 if (A2B <= 0)
2380 goto next;
2381
2382 /* z-coordinate of intersection
2383 */
2384 temp = sqrt(A2B);
2385 t_1 = ((-A) - temp);
2386 t_2 = ((-A) + temp);
2387
2388 t_1_z = t_1 / Tani;
2389 t_2_z = t_2 / Tani;
2390
2391 Zmin = MIN(t_1_z,t_2_z) + z_z;
2392
2393 if (Zmin > HeightPuff)
2394 goto next;
2395
2396 Zmax = MAX(t_1_z,t_2_z) + z_z;
2397
2398 if (Zmax < (-HeightPuff))
2399 goto next;
2400
2401 if (HeightPuff != HeightOut)
2402 {
2403 if (Sini > FLT_EPSILON)
2404 {
2405 double a, hs_dist, t_xy, Z;
2406
2407 Z = t_1_z + z_z;
2408
2409 if ((Z > HeightOut && Z <= HeightPuff) ||
2410 (Z < (-HeightOut) && Z >= (-HeightPuff)))
2411 {
2412 t_xy = t_1 / Sini;
2413 a = toAngle(1.0 - (x_z + t_xy*lzero), -(y_z + t_xy*mzero), Rout);
2414 hs_dist = DIST_CIRCLE(a, hs_long);
2415
2416 if (hs_dist <= hs_ext)
2417 {
2418 if (obscuredByPuff(fabs(hs_dist/hs_ext), Z, HeightPuff))
2419 {
2420 goto set_eclipsed;
2421 }
2422 }
2423 }
2424
2425 Z = t_2_z + z_z;
2426
2427 if ((Z > HeightOut && Z <= HeightPuff) ||
2428 (Z < (-HeightOut) && Z >= (-HeightPuff)))
2429 {
2430 t_xy = t_2 / Sini;
2431 a = toAngle(1.0 - (x_z + t_xy*lzero), -(y_z + t_xy*mzero), Rout);
2432 hs_dist = DIST_CIRCLE(a, hs_long);
2433
2434 if (hs_dist <= hs_ext)
2435 {
2436 if (obscuredByPuff(fabs(hs_dist/hs_ext), t_2_z + z_z, HeightPuff))
2437 {
2438 goto set_eclipsed;
2439 }
2440 }
2441 }
2442 }
2443
2444 if (Zmin > HeightOut)
2445 goto next;
2446
2447 if (Zmax < (-HeightOut))
2448 goto next;
2449 }
2450
2451 /* Compute intersection with inner cylinder
2452 */
2453 B = R0 - (Rin * Rin);
2454 A2B = (A*A) - B;
2455
2456 if (A2B > 0)
2457 {
2458 /* z-coordinate of intersection
2459 */
2460 temp = sqrt(A2B);
2461
2462 t_1 = ((-A) - temp) / Tani;
2463 t_2 = ((-A) + temp) / Tani;
2464
2465 Zmin = MIN(t_1,t_2) + z_z;
2466
2467 if (Zmin < (-HeightIn))
2468 {
2469 Zmax = MAX(t_1,t_2) + z_z;
2470
2471 if (Zmax > HeightIn)
2472 goto next;
2473 }
2474 }
2475
2476 set_eclipsed:
2477 SurfPtr->EclipseFlag = 1;
2478 SurfPtr->visibility = 0;
2479
2480 next:
2481 ++SurfPtr;
2482 ++count;
2483 }
2484
2485 return;
2486 }
2487 #endif
2488