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 <stdlib.h>
22 #include <stdio.h>
23 #include <string.h>
24 #include <ctype.h>
25 #include <float.h>
26 #include "Light.h"
27
28 /* contains miscellaneous small routines */
29
30 /*------------------------------------------------------------------------
31 Copyright (c) 2000, Michael P.D. Bramley.
32
33 Permission is granted to use this code without restriction as long as
34 this copyright notice appears in all source files.
35 */
36
DefineAxis(double * Min,double * Max,double * Inc)37 void DefineAxis(double *Min, double *Max, double *Inc)
38 {
39 /* define local variables... */
40
41 double
42 Test_inc, /* candidate increment value */
43 Test_min, /* minimum scale value */
44 Test_max, /* maximum scale value */
45 Range = *Max - *Min ; /* range of data */
46
47 int i = 0 ; /* counter */
48
49 /* don't create problems -- solve them */
50
51 if( Range < 0 ) {
52 *Inc = 0 ;
53 return ;
54 }
55
56 /* handle special case of repeated values */
57
58 else if( Range == 0) {
59 *Min = ceil(*Max) - 1 ;
60 *Max = *Min + 1 ;
61 *Inc = 1 ;
62 return ;
63 }
64
65 /* compute candidate for increment */
66
67 Test_inc = pow(10, ceil( log10( Range/10 ) ) ) ;
68
69 /* establish maximum scale value... */
70
71 Test_max = ( (long)(*Max / Test_inc)) * Test_inc ;
72
73 if( Test_max < *Max )
74 Test_max += Test_inc ;
75
76 /* establish minimum scale value... */
77
78 Test_min = Test_max ;
79
80 do {
81 ++i ;
82 Test_min -= Test_inc ;
83 }
84 while( Test_min > *Min ) ;
85
86 /* subtracting small values can screw up the scale limits,
87 eg: if DefineAxis is called with (min,max)=(0.01, 0.1),
88 then the calculated scale is 1.0408E17 TO 0.05 BY 0.01.
89 the following if statment corrects for this... */
90
91 if( fabs(Test_min) < 1E-10 )
92 Test_min = 0 ;
93
94 /* adjust for too few tick marks */
95
96 if( i < 6 ) {
97 Test_inc /= 2 ;
98 if( (Test_min + Test_inc) <= *Min )
99 Test_min += Test_inc ;
100 if( (Test_max - Test_inc) >= *Max )
101 Test_max -= Test_inc ;
102 }
103
104 /* pass back axis definition to caller */
105
106 *Min = Test_min ;
107 *Max = Test_max ;
108 *Inc = Test_inc ;
109 }
110
111 /******************************************************************
112 @package nightfall
113 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
114 @version 1.70
115 @short Compute internal hot spot longitude
116 @param (void)
117 @return (double) Internal longitude
118 @heading Disk Setup
119 *******************************************************************/
internal_longitude()120 double internal_longitude()
121 {
122 double ret = -(Binary[Disk].longitudeHS + 180);
123 ret = ((ret < 0) ? (360.0 + ret) : ret);
124 ret = fmod(DTOR * ret, 2.0 * PI);
125
126 return ret;
127 }
128
129 /******************************************************************
130 @package nightfall
131 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
132 @version 1.70
133 @short Map angle in degree into [0, 2*PI)
134 @param (void)
135 @return (double) Radian
136 @heading Disk Setup
137 *******************************************************************/
dtor_2pi(double angle)138 double dtor_2pi(double angle)
139 {
140 double ret = angle;
141
142 ret = ((ret < 0) ? (2.0 * PI + ret) : ret);
143 ret = fmod(DTOR * ret, 2.0 * PI);
144
145 return ret;
146 }
147
148 /******************************************************************
149 @package nightfall
150 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
151 @version 1.0
152 @short Get data path
153 @param (void)
154 @return (char) * The path
155 @heading GUI
156 *******************************************************************/
data_data_fls()157 char * data_data_fls()
158 {
159 static char fls_path[1024];
160 int i;
161
162 for (i=0; i<1024; ++i) fls_path[i] = '\0';
163
164 if (getenv("NIGHTFALL_DATA_DIR") != NULL) {
165 strncpy(fls_path, getenv("NIGHTFALL_DATA_DIR"), 1023);
166 return fls_path;
167 }
168
169 if (getenv("NIGHTFALL_DATAROOT") != NULL) {
170 strncpy(fls_path, getenv("NIGHTFALL_DATAROOT"), (1023 - 5) );
171 strcat(fls_path, "/data");
172 return fls_path;
173 }
174
175 strncpy(fls_path, DATAFLS, (1023 - 5) );
176 return fls_path;
177 }
178
179 /******************************************************************
180 @package nightfall
181 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
182 @version 1.0
183 @short Get cfg path
184 @param (void)
185 @return (char) * The path
186 @heading GUI
187 *******************************************************************/
data_cfg_fls()188 char * data_cfg_fls()
189 {
190 static char fls_path[1024];
191 int i;
192
193 for (i=0; i<1024; ++i) fls_path[i] = '\0';
194
195 if (getenv("NIGHTFALL_CFG_DIR") != NULL) {
196 strncpy(fls_path, getenv("NIGHTFALL_CFG_DIR"), 1023);
197 return fls_path;
198 }
199
200 if (getenv("NIGHTFALL_DATAROOT") != NULL) {
201 strncpy(fls_path, getenv("NIGHTFALL_DATAROOT"), (1023 - 5) );
202 strcat(fls_path, "/cfg");
203 return fls_path;
204 }
205
206 strncpy(fls_path, CFGFLS, (1023 - 5) );
207 return fls_path;
208 }
209
210 /******************************************************************
211 @package nightfall
212 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
213 @version 1.0
214 @short Get doc path
215 @param (void)
216 @return (char) * The path
217 @heading GUI
218 *******************************************************************/
data_doc_fls()219 char * data_doc_fls()
220 {
221 static char fls_path[1024];
222 int i;
223
224 for (i=0; i<1024; ++i) fls_path[i] = '\0';
225
226 if (getenv("NIGHTFALL_DOC_DIR") != NULL) {
227 strncpy(fls_path, getenv("NIGHTFALL_DOC_DIR"), 1023);
228 return fls_path;
229 }
230
231 if (getenv("NIGHTFALL_DATAROOT") != NULL) {
232 strncpy(fls_path, getenv("NIGHTFALL_DATAROOT"), (1023 - 5) );
233 strcat(fls_path, "/doc");
234 return fls_path;
235 }
236
237 strncpy(fls_path, DOCFLS, (1023 - 5) );
238 return fls_path;
239 }
240
241 /******************************************************************
242 @package nightfall
243 @author Markus Kuster (kuster@astro.uni-tuebingen.de)
244 @version 1.0
245 @short Get pixmap path
246 @param (void)
247 @return (char) * The path
248 @heading GUI
249 *******************************************************************/
data_pix_fls()250 char * data_pix_fls()
251 {
252 static char fls_path[1024];
253 int i;
254
255 for (i=0; i<1024; ++i) fls_path[i] = '\0';
256
257 if (getenv("NIGHTFALL_PIX_DIR") != NULL) {
258 strncpy(fls_path, getenv("NIGHTFALL_PIX_DIR"), 1023);
259 return fls_path;
260 }
261
262 if (getenv("NIGHTFALL_DATAROOT") != NULL) {
263 strncpy(fls_path, getenv("NIGHTFALL_DATAROOT"), (1023 - 5) );
264 strcat(fls_path, "/pixmaps");
265 return fls_path;
266 }
267
268 strncpy(fls_path, PIXFLS, (1023 - 5) );
269 return fls_path;
270 }
271
272 static int gd_size = 32;
273
274 static float gd_tt[] = {
275 2012.6, 2322.9, 2646.3, 2824.5, 3054.2,
276 3238.7, 3302.6, 3434.3, 3548.1, 3665.6,
277 3938, 4230.6, 4634.8, 5077.5, 5562.5,
278 5709.5, 6093.9, 6254.9, 6504.3, 6676,
279 7033.3, 7125.6, 7313.8, 7506.9, 7605.4,
280 7705.2, 7806.3, 7908.7, 8012.4, 8117.5,
281 8835.2, 10467, };
282
283 static float gd_beta[] = {
284 0.22016, 0.22326, 0.22016, 0.2031, 0.1938,
285 0.1845, 0.17054, 0.18295, 0.1845, 0.1938,
286 0.23256, 0.34419, 0.4155, 0.43411, 0.4031,
287 0.3876, 0.34574, 0.32248, 0.29922, 0.27287,
288 0.24806, 0.2031, 0.15814, 0.17364, 0.31318,
289 0.58605, 0.72868, 0.8155, 0.93643, 0.94574,
290 0.9876, 1, };
291
292 /******************************************************************
293 @package nightfall
294 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
295 @version 1.0
296 @short Compute gravitational darkening exponent
297 @param (float) Temperature
298 @return (float) Beta, the gravitational darkening exponent
299 @heading Light Curve
300 *******************************************************************/
getGravDarkExp(float temp)301 float getGravDarkExp (float temp)
302 {
303 int i;
304 float val = 1.0;
305
306 gd_beta[gd_size-1] = 1.0;
307
308 if (temp <= gd_tt[0])
309 return gd_beta[0];
310 if (temp >= gd_tt[gd_size-1])
311 return 1.0;
312 for (i = 1; i < gd_size; ++i)
313 {
314 if (fabs(gd_tt[i] - temp) < FLT_EPSILON)
315 return gd_beta[i];
316
317 if (gd_tt[i] > temp)
318 {
319 val = gd_beta[i-1] * ((gd_tt[i]-temp)/(gd_tt[i] - gd_tt[i-1]))
320 + gd_beta[i] * ((temp-gd_tt[i-1])/(gd_tt[i] - gd_tt[i-1]));
321 return val;
322 }
323 }
324 return val;
325 }
326
327
328 static float Albedo[NUM_COMP] = { 0.0 };
329 #define ALBEDO_OFFSET 1.0
330
331 /* The purpose of ALBEDO_OFFSET is to verify whether Albedo[] has
332 * been set manually -facepalms-
333 */
SetAlbedo(int Comp,const char * str)334 void SetAlbedo (int Comp, const char * str)
335 {
336 float ff = atof(str);
337 if (Comp >=0 && Comp < NUM_COMP && 0.0 <= ff && ff <= 1.0)
338 Albedo[Comp] = (ff + ALBEDO_OFFSET);
339 else
340 WARNERR(_("Bad albedo value: "), str);
341 return;
342 }
343
isAlbedoSet(int Comp)344 int isAlbedoSet(int Comp)
345 {
346 if (Albedo[Comp] >= ALBEDO_OFFSET)
347 return 1;
348 return 0;
349 }
350
GetAlbedo(int Comp)351 float GetAlbedo(int Comp)
352 {
353 return (Albedo[Comp] - ALBEDO_OFFSET);
354 }
355
ComputeGravDark()356 void ComputeGravDark()
357 {
358 #ifdef HAVE_DISK
359 int ncomp = 3;
360 #else
361 int ncomp = 2;
362 #endif
363
364 int i;
365 double limit = 7700.0;
366
367 if (getenv("NIGHTFALL_RADIATIVE") != NULL)
368 limit = atof(getenv("NIGHTFALL_RADIATIVE"));
369 if (limit <= 4000. || limit >= 12000.) limit = 7700.0;
370
371 /* >>>>>>>>>>> determine albedo and GravDarkCoeff <<<<<<<< */
372
373 for (i = 0; i < ncomp; ++i) {
374
375 if (Binary[i].Temperature <= limit) {
376
377 /* -------------- convective ---------------------- */
378 if (Albedo[i] >= ALBEDO_OFFSET)
379 Binary[i].Albedo = Albedo[i] - ALBEDO_OFFSET;
380 else
381 Binary[i].Albedo = 0.5;
382
383 Binary[i].GravDarkCoeff = getGravDarkExp(Binary[i].Temperature);
384 Binary[i].GravDarkCoeff = Binary[i].GravDarkCoeff / 4.0;
385
386 } else {
387
388 /* ---------------- radiative ------------------------ */
389 if (Albedo[i] >= ALBEDO_OFFSET)
390 Binary[i].Albedo = Albedo[i] - ALBEDO_OFFSET;
391 else
392 Binary[i].Albedo = 1.0;
393
394 Binary[i].GravDarkCoeff = 0.25;
395 }
396
397 if (Flags.debug[VERBOSE] == ON)
398 printf(_(" Comp %2d GravDarkCoeff %8.4f Albedo %8.4f\n"),
399 i+1, Binary[i].GravDarkCoeff, Binary[i].Albedo);
400 }
401
402 return;
403
404 }
405
406 /****************************************************************************
407 @package nightfall
408 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
409 @version 1.0
410 @short Set Xwindow size for PGPLOT
411 @tip Allocates 21 * sizeof(char) without freeing (called only once)
412 @param (void)
413 @return (void)
414 @heading Plotting
415 ****************************************************************************/
SetPgplotWindow()416 void SetPgplotWindow()
417 {
418 int width = 551; /* width and height */
419 char *windowsize; /* the window size */
420 float frac; /* fractional display width */
421 char * frac_size; /* fractional display width */
422 static char hasbuf[32];
423 char gnu_geometry[256];
424
425 /* this function is only called at startup, we never free frac_size */
426 /* unless program exits (not really memory leak) */
427
428 /*@i@*/ frac_size = malloc(21);
429 if (frac_size == NULL) nf_error(_(errmsg[0]));
430
431 /* check whether environment variable is defined */
432
433 windowsize = getenv("PGPLOT_XW_WIDTH");
434
435 /* if not, set fractional display size */
436 if (windowsize == NULL) {
437
438 if (getenv("GNUPLOT_GEOMETRY") == NULL)
439 strncpy(gnu_geometry, GNU_GEOMETRY, 255);
440 else
441 strncpy(gnu_geometry, getenv("GNUPLOT_GEOMETRY"), 255);
442 gnu_geometry[255] = '\0';
443
444 (void) sscanf(gnu_geometry, "%dx", &width);
445 frac = (float)width/867.0; /* 867 pix is builtin default for PGPLOT */
446 frac = frac / 1.2;
447 sprintf(frac_size,"PGPLOT_XW_WIDTH=%4.2f", frac);
448 #ifdef HAVE_PUTENV
449 /* is putenv() POSIX ? */
450 /*@i@*/ putenv(frac_size);
451 #else
452 free(frac_size);
453 #endif
454 }
455
456 if (NULL == getenv("PGPLOT_BUFFER")) {
457 sprintf(hasbuf,"PGPLOT_BUFFER=yes");
458 putenv(hasbuf);
459 }
460 /*@i@*/ return;
461
462 }
463
464 /****************************************************************************
465 @package nightfall
466 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
467 @version 1.0
468 @short Print error message and exit
469 @param (char) error[] The error text
470 @return (void)
471 @heading Error handling (non-interactive)
472 ****************************************************************************/
nf_error(char * error_msg)473 void nf_error(char *error_msg)
474 {
475 int i = 0, status = 0;
476
477 char command[256+64];
478 char command_one[256] = " \n";
479
480 while (i < 255 && *error_msg != '\0') {
481 command_one[i] = *error_msg;
482 ++error_msg;
483 ++i;
484 }
485
486 status =
487 system("which xmessage > /dev/null");
488
489 if (status == 0) {
490 sprintf(command,
491 _("xmessage \"ERROR: %s\" "), command_one);
492 status =
493 system(command);
494 }
495
496 fprintf(stderr, _("**ERROR**: %s\n"), command_one);
497 fprintf(stderr, _(".. exit to system\n"));
498
499 exit (EXIT_FAILURE);
500 }
501
502
503 /****************************************************************************
504 @package nightfall
505 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
506 @version 1.0
507 @short Convert string to lowercase
508 @tip strlwr is not ANSI
509 @param (char) *instring The input string
510 @return (void)
511 @heading String Handling
512 ****************************************************************************/
nf_strlwr(char * instring)513 void nf_strlwr(char * instring)
514 {
515 if (instring == NULL) return;
516
517 while ( (*instring) != '\0') {
518 if (isalpha(*instring) != 0) (*instring) = tolower(*instring);
519 ++instring;
520 }
521 }
522
523
524 /****************************************************************************
525 @package nightfall
526 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
527 @version 1.0
528 @short Allocate a double 2D matrix
529 @tip Code snippet from comp.lang.c FAQ, modified
530 @param (long) nrh number of rows
531 @param (long) nch number of collumns
532 @return (double**) pointer to allocated memory
533 @heading Memory Handling
534 ****************************************************************************/
dmatrix(long nrows,long ncols)535 double **dmatrix(long nrows, long ncols)
536 {
537 long i;
538 double **m;
539
540 /* allocate pointers to rows */
541 /*@ignore@*/
542 /* cast argument of malloc b/o NEC SX-5 compiler warning */
543 m = malloc((size_t) (nrows) * sizeof(double*));
544 if (m == NULL) return (NULL);
545
546 /* allocate rows and set pointers to them */
547 for (i = 0; i < nrows; ++i) {
548 m[i] = malloc((size_t)(ncols)*sizeof(double));
549 if (!m[i]) return (NULL);
550 }
551 /* return pointer to array of pointers to rows */
552 return m;
553 /*@end@*/
554 }
555
556 /****************************************************************************
557 @package nightfall
558 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
559 @version 1.0
560 @short Free a double 2D matrix
561 @param (double**) 2D matrix allocated with dmatrix()
562 @param (long) nrh number of rows
563 @param (long) nch number of collumns
564 @return (void)
565 @heading Memory Handling
566 ****************************************************************************/
free_dmatrix(double ** m,long nrows,long ncols)567 void free_dmatrix(double **m, long nrows, long ncols)
568 {
569 long i;
570 (void) ncols;
571
572 for (i = 0; i < nrows; ++i) {
573 free(m[i]);
574 }
575
576 free(m);
577 return;
578 }
579
580
581 /****************************************************************************
582 @package nightfall
583 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
584 @version 1.0
585 @short Set normalization phase for output flux (from data file)
586 @param (int) band the bandpass
587 @param (int) phase the phase at which to normalize
588 @return (void)
589 @heading Light Curve
590 ****************************************************************************/
591 static float NormalPhase[NUM_MAG] = {
592 -1.0, -1.0, -1.0, -1.0,
593 -1.0, -1.0, -1.0, -1.0,
594 -1.0, -1.0, -1.0, -1.0
595 };
596
597 double P90[NUM_MAG] = {
598 0.75, 0.75, 0.75, 0.75,
599 0.75, 0.75, 0.75, 0.75,
600 0.75, 0.75, 0.75, 0.75
601 };
602
SetNormalPhase(int band,float phase)603 void SetNormalPhase (int band, float phase)
604 {
605 /* internal consistency check
606 */
607 if (band >= NUM_MAG || band < 0)
608 {
609 if (Flags.debug[VERBOSE] == ON || Flags.debug[WARN] == ON)
610 fprintf(stderr,
611 _("**Warning**: NormalPhase(): Bandpass (%d) out of range (0, %d)\n"),
612 band, NUM_MAG-1);
613 return;
614 }
615
616 /* map phase to 0..1
617 *
618 * if (phase > 1.0)
619 * phase = phase - ((int) phase);
620 * else if (phase < -1.0)
621 * phase = phase - ((int) phase);
622 * phase = (phase < 0.0) ? (phase + 1.0) : phase;
623 */
624
625 /* No negative phases
626 */
627 if (phase < -1.0)
628 phase = phase - ((int) phase);
629 phase = (phase < 0.0) ? (phase + 1.0) : phase;
630
631 NormalPhase[band] = phase;
632
633 return;
634 }
635
636 /****************************************************************************
637 @package nightfall
638 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
639 @version 1.0
640 @short Reset normalization phase for output flux
641 @param (void)
642 @return (void)
643 @heading Light Curve
644 ****************************************************************************/
ResetNormalPhase()645 void ResetNormalPhase ()
646 {
647 register int i;
648
649 for (i = 0; i < NUM_MAG; ++i)
650 NormalPhase[i] = -1.0;
651 return;
652 }
653
ResetNormalPhaseOne(int band)654 void ResetNormalPhaseOne (int band)
655 {
656 if (band >= 0 && band < NUM_MAG)
657 NormalPhase[band] = -1.0;
658 return;
659 }
660
661 /****************************************************************************
662 @package nightfall
663 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
664 @version 1.0
665 @short Set normalization points for output flux
666 @param (int) band the bandpass
667 @param (int) phase the phase at which to normalize
668 @return (void)
669 @heading Light Curve
670 ****************************************************************************/
671 static int NormalPoint[NUM_MAG] = { 0 };
672
NormalPoints()673 void NormalPoints ()
674 {
675 register int j, i;
676 float mindist, dist;
677 float phase;
678
679 for (i = 0; i < NUM_MAG; ++i)
680 {
681 /* NormalPhase is the user-supplied (in data file) phase
682 * at which normalization should occur. Ideally this should be
683 * at quadrature.
684 */
685 if (NormalPhase[i] < 0.0)
686 {
687 NormalPoint[i] = 0;
688 phase = (FluxOut[0].Phase / (PI+PI)) - 0.5;
689
690 /* map phase to 0..1
691 */
692 if (phase > 1.0)
693 phase = phase - ((int) phase);
694 else if (phase < -1.0)
695 phase = phase - ((int) phase);
696 phase = (phase < 0.0) ? (phase + 1.0) : phase;
697
698 P90[i] = phase;
699
700 }
701 else
702 {
703 mindist = NumPhases;
704
705 for (j = 0; j < N_FluxOutFull; ++j)
706 {
707 phase = (FluxOut[j].Phase / (PI+PI)) - 0.5;
708
709 if (NumPhases == 1)
710 {
711 /* map phase to 0..1
712 */
713 if (phase > 1.0)
714 phase = phase - ((int) phase);
715 else if (phase < -1.0)
716 phase = phase - ((int) phase);
717 phase = (phase < 0.0) ? (phase + 1.0) : phase;
718 }
719 else
720 {
721 /* No negative phases
722 */
723 if (phase < -1.0)
724 phase = phase - ((int) phase);
725 phase = (phase < 0.0) ? (phase + 1.0) : phase;
726 }
727
728 /* search for better minimum
729 */
730 if (mindist > (dist = fabs(phase - NormalPhase[i])))
731 {
732 mindist = dist;
733 NormalPoint[i] = j;
734 P90[i] = phase;
735 }
736 }
737 }
738 }
739 return;
740 }
741
742 /****************************************************************************
743 @package nightfall
744 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
745 @version 1.37
746 @short Luminosity for components
747 @param (int) component the bandpass
748 @param (int) phaseindex the current phase index
749 @return (void)
750 @heading Light Curve
751 ****************************************************************************/
LightLuminosity(int component,int phaseindex,int n_phase)752 void LightLuminosity(int component, int phaseindex, int n_phase)
753 {
754 static int check = -1;
755 int band;
756
757 /* Initialize
758 */
759 if (component < 0)
760 {
761 int i;
762
763 check = -1;
764
765 for (i = 0; i < NUM_COMP; ++i) {
766 for (band = 0; band < NUM_MAG; ++band) {
767 D90[i][band] = NumPhases;
768 L90[i][band] = 1.0;
769 }
770 }
771 return;
772 }
773
774 /* check whether all normalization points are -1 */
775 if (check == -1)
776 {
777 check = 0;
778
779 for (band = 0; band < NUM_MAG; ++band)
780 {
781 D90[component][band] = NumPhases;
782
783 if ((NormalPhase[band] != -1))
784 {
785 check = 1;
786 }
787 }
788 }
789
790 if ((check == 0) && (phaseindex == 0))
791 {
792 /* Normalization point is phase index 0 for all bands
793 */
794 for (band = 0; band < NUM_MAG; ++band)
795 {
796 L90[component][band] = FluxOut[0].Mag[band];
797 if (component == Secondary)
798 L90[component][band] -= L90[Primary][band];
799 #ifdef HAVE_DISK
800 if (component == Disk)
801 L90[component][band] -= (L90[Primary][band]+L90[Secondary][band]);
802 #endif
803 }
804 }
805
806 if (check == 1)
807 {
808 float mindist, dist;
809 float phase;
810 int j, jnorm = -1;
811
812 for (band = 0; band < NUM_MAG; ++band)
813 {
814 mindist = D90[component][band];
815
816 /* find best normalization phase jnorm
817 */
818 if (NormalPhase[band] != -1.0)
819 {
820 for (j = 0; j <= phaseindex; ++j)
821 {
822 phase = (FluxOut[j].Phase / (PI+PI)) - 0.5;
823
824 /* map phase to 0..1
825 */
826 if (phase > 1.0)
827 phase = phase - ((int) phase);
828 else if (phase < -1.0)
829 phase = phase - ((int) phase);
830 phase = (phase < 0.0) ? (phase + 1.0) : phase;
831
832 /* Add current phase (orbit) count if we compute
833 * multiple orbits
834 */
835 phase += n_phase;
836
837 /* search for better minimum
838 */
839 if (mindist > (dist = fabs(phase - NormalPhase[band])))
840 {
841 mindist = dist;
842 jnorm = j;
843 }
844 }
845 }
846 else
847 {
848 jnorm = 0;
849 }
850
851 /* copy flux to L90 if we are at best phase
852 */
853 if (jnorm == phaseindex)
854 {
855 D90[component][band] = mindist;
856
857 L90[component][band] = FluxOut[jnorm].Mag[band];
858 if (component == Secondary)
859 L90[component][band] -=
860 L90[Primary][band];
861 #ifdef HAVE_DISK
862 if (component == Disk)
863 L90[component][band] -=
864 (L90[Primary][band]+L90[Secondary][band]);
865 #endif
866 }
867 }
868 }
869
870 return;
871 }
872
873 /****************************************************************************
874 @package nightfall
875 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
876 @version 1.0
877 @short Normalize output flux
878 @param (void)
879 @return (void)
880 @heading Light Curve
881 ****************************************************************************/
NormalizeFlux()882 void NormalizeFlux()
883 {
884 int j, band;
885 double ThirdAdd[NUM_MAG];
886
887 /* ---------------- initialize -------------------------------------- */
888
889 NormalPoints ();
890
891 for (band = 0; band < NUM_MAG; ++band)
892 {
893 ThirdAdd[band] = 0.0;
894 F90[band] = FluxOutFull[NormalPoint[band]].Mag[band];
895 if ( ( 1.0 - Orbit.Third[band]) >= FLT_EPSILON )
896 {
897 ThirdAdd[band] = F90[band]
898 * ( Orbit.Third[band] / ( 1.0 - Orbit.Third[band]) );
899 F90[band] = F90[band] + ThirdAdd[band];
900 }
901 }
902
903 /* ---------------- Normalize Output ------------------------------------ */
904
905 for (j = 0; j < N_FluxOutFull; ++j)
906 {
907 for (band = 0; band < NUM_MAG; ++band)
908 {
909 FluxOutFull[j].Mag[band] = -2.5
910 * log10((ThirdAdd[band] + FluxOutFull[j].Mag[band])/F90[band]);
911 }
912 }
913
914 return;
915 }
916
917 /****************************************************************************
918 @package nightfall
919 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
920 @version 1.0
921 @short Shift phase to correct interval
922 @param (void)
923 @return (void)
924 @heading Light Curve
925 ****************************************************************************/
PhaseShift()926 void PhaseShift()
927 {
928 int i, test;
929
930 test = ON;
931
932 do {
933 if ( (FluxOut[0].Phase - PI) >= FLT_EPSILON ) {
934 for (i = 0; i < PhaseSteps; ++i) {
935 FluxOut[i].Phase=FluxOut[i].Phase-2.0*PI;
936 }
937 test = ON; /* there was something to do */
938 } else {
939 if ( (FluxOut[0].Phase - (-PI)) <= (-FLT_EPSILON) ) {
940 for (i = 0; i < PhaseSteps; ++i) {
941 FluxOut[i].Phase=FluxOut[i].Phase+2.0*PI;
942 }
943 test = ON; /* there was something to do */
944 } else { test = OFF; }
945 }
946 } while (test == ON);
947
948 return;
949 }
950
951 /****************************************************************************
952 @package nightfall
953 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
954 @version 1.0
955 @short Read a single line from an open file
956 @param (char*) s The buffer string
957 @param (int) lim Size of buffer string
958 @param (FILE*) fpe Input filehandle
959 @return (int) Length of line read
960 @heading Light Curve
961 ****************************************************************************/
LireLigne(char * s,int lim,FILE * fpe)962 int LireLigne(char *s, int lim, FILE *fpe)
963 {
964 if (fgets(s, lim, fpe) == NULL){
965 return 0;
966 }
967
968 return (int) strlen(s);
969
970 }
971
972 /****************************************************************************
973 @package nightfall
974 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
975 @version 1.0
976 @short Reset the data file list
977 @param (void)
978 @return (void)
979 @heading Data Fitting
980 ****************************************************************************/
ClearList()981 void ClearList()
982 {
983 FileType *item_ptr;
984 FileType *current;
985
986 item_ptr = FileList;
987 current = FileList;
988
989 while (current != NULL) {
990 item_ptr = current->nextFile;
991 free(current);
992 current = item_ptr;
993 }
994 FileList = NULL;
995
996 /* Reset all normalization points
997 */
998 ResetNormalPhase ();
999
1000 return;
1001 }
1002
1003
1004 /****************************************************************************
1005 @package nightfall
1006 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1007 @version 1.0
1008 @short Add an item to the data file list
1009 @param (char*) s The item
1010 @return (void)
1011 @heading Data Fitting
1012 ****************************************************************************/
AddToList(const char * item,int Format)1013 void AddToList(const char * item, int Format)
1014 {
1015 FileType *new_ptr;
1016
1017 /*@ignore@*/
1018 new_ptr = malloc(sizeof(FileType));
1019
1020 if (new_ptr != NULL) {
1021 strncpy(new_ptr->DataFile, item, MAX_INFILE_PATH);
1022 new_ptr->DataFormat = Format;
1023 new_ptr->nextFile = FileList;
1024 FileList = new_ptr;
1025 }
1026
1027 return;
1028 /*@end@*/
1029 }
1030
1031
1032 /****************************************************************************
1033 @package nightfall
1034 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1035 @version 1.0
1036 @short Potential in Y direction at LagrangeOne
1037 @tip Uses Global Variable Mq (Mass Ratio)
1038 @param (double) y Distance in y-direction
1039 @return (double)
1040 @heading Light Curve
1041 ****************************************************************************/
RocheYPerpendL1(double y)1042 double RocheYPerpendL1(double y)
1043 {
1044 /* Uses Global Variable Mq (Mass Ratio) */
1045 /* Uses Global Variable RochePot */
1046 /* Potential in Y direction at LagrangeOne */
1047 /* used to define the throat */
1048
1049 double Eval; /* return value */
1050 double r2, r1; /* radii from Prim./Second. */
1051 /* pointer */
1052 BinaryComponent *Bptr = &Binary[Primary];
1053
1054 r1 = sqrt(SQR(Bptr->RLag1) + (y*y));
1055 r2 = sqrt(SQR(1.0-Bptr->RLag1) + (y*y));
1056 Eval = -RochePot
1057 + (1.0/r1)
1058 + (Mq * (1.0/r2 - Bptr->RLag1))
1059 + ((Mq+1.0)/2.0) * (SQR(Bptr->RLag1) + (y*y));
1060 return(Eval);
1061 }
1062
1063 /****************************************************************************
1064 @package nightfall
1065 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1066 @version 1.0
1067 @short Potential in Z direction at LagrangeOne
1068 @tip Uses Global Variable Mq (Mass Ratio)
1069 @param (double) z Distance in z-direction
1070 @return (double)
1071 @heading Light Curve
1072 ****************************************************************************/
RocheZPerpendL1(double z)1073 double RocheZPerpendL1(double z)
1074 {
1075 /* Uses Global Variable Mq (Mass Ratio) */
1076 /* Uses Global Variable RochePot */
1077 /* Potential in Z direction at LagrangeOne */
1078 /* used to define the throat */
1079
1080 double Eval; /* return value */
1081 double r2, r1; /* radii from Prim./Second. */
1082 /* pointer */
1083 BinaryComponent *Bptr = &Binary[Primary];
1084
1085 r1 = sqrt( (Bptr->RLag1*Bptr->RLag1) + (z*z));
1086 r2 = sqrt( (1.0-Bptr->RLag1)*(1.0-Bptr->RLag1) + (z*z));
1087 Eval = -RochePot
1088 + (1.0/r1)
1089 + (Mq * (1.0/r2 - Bptr->RLag1))
1090 + ((Mq+1.0)/2.0) * (Bptr->RLag1 * Bptr->RLag1 ) ;
1091 return(Eval);
1092 }
1093
1094 /****************************************************************************
1095 @package nightfall
1096 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1097 @version 1.0
1098 @short Potential in Z direction at Secondary
1099 @tip Uses Global Variable Mq (Mass Ratio)
1100 @param (double) z Distance in z-direction
1101 @return (double)
1102 @heading Light Curve
1103 ****************************************************************************/
RochePerpendSecond(double z)1104 double RochePerpendSecond(double z)
1105 {
1106 /* Uses Global Variable Mq (Mass Ratio) */
1107 /* Uses Global Variable RochePot */
1108
1109 double Eval; /* return value */
1110 double r2, r1; /* radii from Prim./Second. */
1111
1112 r1 = sqrt(1.0 + (z*z));
1113 r2 = z;
1114
1115 if (z > DBL_EPSILON)
1116 Eval = -RochePot
1117 + (1.0/r1)
1118 + (Mq * (1.0/r2 - 1.0))
1119 + ((Mq+1.0)/2.0);
1120 else
1121 {
1122 if (Flags.debug[VERBOSE] == ON || Flags.debug[WARN] == ON)
1123 fprintf(stderr, "**Warning**: Divide by zero, File %s, Line %d \n",
1124 __FILE__, __LINE__);
1125 Eval = -RochePot
1126 + (1.0/r1)
1127 + ((Mq+1.0)/2.0);
1128 }
1129 return(Eval);
1130 }
1131
1132 /****************************************************************************
1133 @package nightfall
1134 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1135 @version 1.0
1136 @short Potential
1137 @tip Uses Global Variables Mq, RochePot, lambda, nu, F
1138 @param (double) x Radius from origin
1139 @return (double)
1140 @heading Light Curve
1141 ****************************************************************************/
RocheSurface(double x)1142 double RocheSurface(double x)
1143 {
1144 /* Uses Global Variable Mq (Mass Ratio) */
1145 /* Uses Global Variable RochePot */
1146 /* Uses Global Variable lambda (Coordinate) */
1147 /* Uses Global Variable nu (Coordinate) */
1148 /* Uses Global Variable F (nonsync ) */
1149
1150 /* input x is radius from origin */
1151
1152 double Eval; /* return value */
1153 double sqrX; /* temp variable */
1154 double lamX; /* temp variable */
1155
1156 sqrX = x * x;
1157 lamX = x*lambda;
1158
1159 Eval = -RochePot
1160 + (1.0/x)
1161 + (Mq * (1.0/sqrt(1.0 + sqrX - 2*lamX) - lamX))
1162 + F*F* ((Mq+1.0)/2.0) * (sqrX * (1.0 - (nu * nu)));
1163 return(Eval);
1164 }
1165
1166 /****************************************************************************
1167 @package nightfall
1168 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1169 @version 1.0
1170 @short Potential in x-y plane
1171 @param (double) x x-vector
1172 @param (double) y y-vector
1173 @param (double) q mass ratio
1174 @param (double) f async rotation
1175 @return (double)
1176 @heading Light Curve
1177 ****************************************************************************/
RocheXYPrimary(double x,double y,double q,double f)1178 double RocheXYPrimary(double x, double y, double q, double f)
1179 {
1180 double Eval; /* return value */
1181 double r1, r2; /* radii from Prim./Second. */
1182 double sqrX; /* temp variable */
1183
1184
1185 sqrX = x * x;
1186 r1 = sqrt( sqrX +(y*y) );
1187 r2 = sqrt( SQR(x-1.0) + (y*y) );
1188
1189 if (r1 > DBL_EPSILON && r2 > DBL_EPSILON)
1190 Eval = 1.0/r1
1191 + q * ( 1.0/r2 - x )
1192 + ( (q+1.0)/2.0 ) * ( sqrX + (y*y) ) * (f*f);
1193 else
1194 {
1195 if (Flags.debug[VERBOSE] == ON || Flags.debug[WARN] == ON)
1196 fprintf(stderr, "**Warning**: Divide by zero, File %s, Line %d \n",
1197 __FILE__, __LINE__);
1198 Eval = 1.0;
1199 }
1200 return(Eval);
1201 }
1202
1203 /****************************************************************************
1204 @package nightfall
1205 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1206 @version 1.0
1207 @short Analytic polar radius
1208 @param (double) q Mass ratio
1209 @param (double) r Mean radius
1210 @return (double)
1211 @heading Light Curve
1212 ****************************************************************************/
PolRad(double q,double r)1213 double PolRad(double q, double r)
1214 {
1215 /* q is mass ratio, r mean radius */
1216
1217 /* these are the legendre functions for argument zero */
1218 static double Pl2 = -1.0/2.0;
1219 static double Pl4 = 3.0/8.0;
1220 static double Pl6 = -5.0/16.0;
1221 static double Pl8 = 35.0/128.0;
1222 static double Pl10 = -315.0/1280.0;
1223
1224 register double powr; /* power of radius */
1225 register double re; /* return value */
1226
1227 powr = r*r*r; /* r^3 */
1228 re = powr*q*Pl2;
1229
1230 powr = powr*r*r; /* r^5 */
1231 re = re + powr*q*Pl4;
1232
1233 powr = powr*r; /* r^6 */
1234 re = re + powr*3.0*q*Pl2*q*Pl2;
1235
1236 powr = powr*r; /* r^7 */
1237 re = re + powr*q*Pl6;
1238
1239 powr = powr*r; /* r^8 */
1240 re = re + powr*8.0*q*q*Pl2*Pl4;
1241
1242 powr = powr*r; /* r^9 */
1243 re = re + powr*q*Pl8 + powr*12.0*(q*Pl2)*(q*Pl2)*(q*Pl2);
1244
1245 powr = powr*r; /* r^10*/
1246 re = re + powr*10.0*q*q*Pl2*Pl6 + powr*5.0*q*q*Pl4*Pl4;
1247
1248 powr = powr*r; /* r^11*/
1249 re = re + powr*q*Pl10 + powr*55.0*q*q*q*Pl2*Pl2*Pl4;
1250
1251 re = r*re +r;
1252
1253 return(re);
1254 }
1255
1256 /****************************************************************************
1257 @package nightfall
1258 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1259 @version 1.0
1260 @short Returns (volume-scaledvolume) as function of r_0
1261 @tip Zero if r_0 = RadMean(ScaledVolume)
1262 @param (double) r_0 Radius
1263 @return (double)
1264 @heading Light Curve
1265 ****************************************************************************/
VolOfR(double r_0)1266 double VolOfR(double r_0)
1267 {
1268 /* uses global variable PhaseScaledVolume */
1269 /* uses global variable Mq */
1270 /* Volume analytic as by Kopal exact to O(r**11) */
1271
1272 register double NN; /* some temp variable */
1273 register double RadPower; /* power of radius */
1274 register double Vol1 = 1.0; /* volume */
1275 double Volume; /* return value */
1276
1277 NN = (Mq + 1.0) / 2.0;
1278
1279 RadPower = r_0 * r_0 * r_0; /* r^3 */
1280 Vol1 = Vol1 + 2.0 * NN * RadPower;
1281
1282 RadPower = SQR(RadPower); /* r^6 */
1283 Vol1 = Vol1 + (12.0/5.0)*(Mq*Mq)*RadPower
1284 + (32.0/5.0)*(NN*NN)*RadPower
1285 + (8.0/5.0)*NN*Mq*RadPower;
1286
1287 RadPower = RadPower*SQR(r_0); /* r^8 */
1288 Vol1 = Vol1 + (15.0/7.0)*(Mq*Mq)*RadPower;
1289
1290 RadPower = RadPower*r_0; /* r^9 */
1291 Vol1 = Vol1 + (22.0/7.0)*(Mq*Mq)*Mq*RadPower
1292 + (176.0/7.0)*(NN*NN)*NN*RadPower
1293 + (296.0/35.0)*NN*Mq
1294 *(2.0*Mq + NN)*RadPower;
1295
1296 RadPower = RadPower*r_0; /* r^10*/
1297 Vol1 = Vol1 + (18.0/9.0)*(Mq*Mq)*RadPower;
1298
1299 RadPower = RadPower*r_0; /* r^11*/
1300 Vol1 = Vol1 + (157.0/7.0)*(Mq*Mq)*Mq*RadPower
1301 + (26.0/35.0)*NN*Mq
1302 *(Mq + 3.0*NN)*RadPower;
1303
1304 Volume = Vol1
1305 * (4.0*PI/3.0)*(r_0*r_0*r_0);
1306
1307 Volume = Volume - PhaseScaledVolume;
1308
1309 return(Volume);
1310
1311 }
1312
1313 /****************************************************************************
1314 @package nightfall
1315 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1316 @version 1.0
1317 @short Zero at LagrangeOne
1318 @param (double) x Position on x axix
1319 @return (double)
1320 @heading Light Curve
1321 ****************************************************************************/
LagrangeOne(double x)1322 double LagrangeOne(double x)
1323 {
1324 /* Uses Global Variable Mq (Mass Ratio) */
1325 /* Eval is Zero at LagrangeOne (F=1) */
1326 /* or at critical Radius (F != 1) */
1327
1328 double Eval = - 1.0/(x*x) + Mq * ( 1.0/SQR(x-1) - 1.0 ) + F*F*x*( Mq + 1);
1329
1330 return(Eval);
1331 }
1332
1333 /****************************************************************************
1334 @package nightfall
1335 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1336 @version 1.0
1337 @short Zero at LagrangeTwo
1338 @param (double) x Position on x axix
1339 @return (double)
1340 @heading Light Curve
1341 ****************************************************************************/
LagrangeTwo(double x)1342 double LagrangeTwo(double x)
1343 {
1344 /* Uses Global Variable Mq (Mass Ratio) */
1345 /* Eval is Zero at LagrangeTwo (F=1) */
1346 /* or at critical Radius (F != 1) */
1347
1348 double Eval = - 1.0/(x*x)
1349 + Mq * ( (1.0-x)/((x-1.0)*SQR(x-1.0)) - 1.0 )
1350 + F*F*x*( Mq + 1);
1351 return(Eval);
1352 }
1353
1354
1355 /****************************************************************************
1356 @package nightfall
1357 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1358 @version 1.0
1359 @short Zero at LagrangeThree
1360 @param (double) x Position on x axix
1361 @return (double)
1362 @heading Light Curve
1363 ****************************************************************************/
LagrangeThree(double x)1364 double LagrangeThree(double x)
1365 {
1366 /* Uses Global Variable Mq (Mass Ratio) */
1367 /* Eval is Zero at LagrangeThree (F=1) */
1368 /* or at critical Radius (F != 1) */
1369
1370 double Eval = - x/((-x)*(-x)*(-x))
1371 + Mq * ( (1.0)/SQR(1.0-x) - 1.0 )
1372 + F*F*x*( Mq + 1);
1373 return(Eval);
1374 }
1375
1376 /****************************************************************************
1377 @package nightfall
1378 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1379 @version 1.0
1380 @short Find root of a function in [Low,Up] with Tolerance
1381 @long Uses bisection algorithm. f(Low) and f(Up) must have opposite signs,
1382 thus indicating that they bracket a root of f().
1383 If this is not the case, the program exits with an error message.
1384 @param (double)(*Function)(double) Pointer to the function
1385 @param (const double) Low Lower limit
1386 @param (const double) Up Upper limit
1387 @param (const double) Tolerance The convergence criterion
1388 @param (const char*) caller The calling function
1389 @param (int*) Err The exit status
1390 @return (double)
1391 @heading Light Curve
1392 ****************************************************************************/
RootBisect(double (* Function)(double),const double Low,const double Up,const double Tolerance,const char * caller,int * Err)1393 double RootBisect(double (*Function)(double),
1394 const double Low, const double Up,
1395 const double Tolerance, const char *caller, int *Err)
1396 {
1397 double a = Low, b = Up, c = Up; /* storage */
1398 double sa, sc; /* sign */
1399 double fa = (*Function)(a); /* function value */
1400 double fb = (*Function)(b); /* function value */
1401 double fc; /* function value */
1402 double Tolerance1;
1403 unsigned int i = 0;
1404
1405 if ( (fa >= FLT_EPSILON && fb >= FLT_EPSILON)
1406 || (fa <= (-FLT_EPSILON) && fb <= (-FLT_EPSILON) )) {
1407
1408 if (Flags.interactive == OFF) {
1409 fprintf(stderr,
1410 _("** ERROR **: %s: RootFind: Root not bracketed \n"), caller);
1411 fprintf(stderr,
1412 _("This is likely due to an extreme combination of\n"));
1413 fprintf(stderr, _("mass ratio and fill factor\n"));
1414 abort();
1415 } else {
1416 fprintf(stderr,
1417 _("** ERROR **: %s: RootFind: Root not bracketed \n"), caller);
1418 fprintf(stderr,
1419 _("This is likely due to an extreme combination of\n"));
1420 fprintf(stderr, _("mass ratio and fill factor\n"));
1421 *Err = 1;
1422 return(0.0);
1423 }
1424 }
1425
1426 do {
1427 c = (a+b)*0.5;
1428
1429 fc = (*Function)(c);
1430
1431 Tolerance1 = (2.0 * FLT_EPSILON * fabs(c)) + (0.5 * Tolerance);
1432
1433 if ( (fabs(a-c) < Tolerance1) || (fabs(fc) < Tolerance) )
1434 return c;
1435
1436 /*
1437 sa = (fa < 0) ? -1.0 : 1.0;
1438 sc = (fc < 0) ? -1.0 : 1.0;
1439 */
1440 sa = (fa >= 0) - (fa < 0);
1441 sc = (fc >= 0) - (fc < 0);
1442 if (sa != sc) { b = c; fb = fc; }
1443 else { a = c; fa = fc; }
1444
1445 ++i;
1446 } while (i < 1000);
1447
1448 return c;
1449 }
1450
1451 /****************************************************************************
1452 @package nightfall
1453 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
1454 @version 1.0
1455 @short Find root of a function in [Low,Up] with Tolerance
1456 @long Uses Brents algorithm. f(Low) and f(Up) must have opposite signs,
1457 thus indicating that they bracket a root of f().
1458 If this is not the case, the program exits with an error message.
1459 Adopted from procedure 'zero' in: Richard Brent, Algorithms for
1460 minimization without derivatives, Prentice-Hall, Inc. (1973)
1461 @param (double)(*Function)(double) Pointer to the function
1462 @param (const double) Low Lower limit
1463 @param (const double) Up Upper limit
1464 @param (const double) Tolerance The convergence criterion
1465 @param (const char*) caller The calling function
1466 @param (int*) Err The exit status
1467 @return (double)
1468 @heading Light Curve
1469 ****************************************************************************/
RootFind(double (* Function)(double),const double Low,const double Up,const double Tolerance,const char * caller,int * Err)1470 double RootFind(double (*Function)(double),
1471 const double Low, const double Up,
1472 const double Tolerance, const char *caller, int *Err)
1473 {
1474
1475 register unsigned long i = 0; /* loop variable */
1476 double a = Low, b = Up, c = Up; /* storage */
1477 double Delta = 0.0, e = 0.0; /* convergence testing */
1478 double test1, test2; /* convergence testing */
1479 double fa = (*Function)(a); /* function value */
1480 double fb = (*Function)(b); /* function value */
1481 double fc; /* function value */
1482 double p, q, r, s; /* for quadratic interpol. */
1483 double Tolerance1 = 0.0; /* convergence testing */
1484 double Bisect, Eval; /* bisection */
1485 static double laststore = 0.5; /* last root found */
1486 register unsigned long j = 1; /* loop variable */
1487 int is_ok = 0; /* brackets found */
1488
1489 *Err = 0;
1490
1491
1492 /* -------- whether interval (a,b) brackets a root ------------------- */
1493
1494 if ( (fa >= FLT_EPSILON && fb >= FLT_EPSILON)
1495 || (fa <= (-FLT_EPSILON) && fb <= (-FLT_EPSILON) )) {
1496
1497 /* root not bracketed, try last value */
1498
1499 WARNING(_("Root not bracketed, try to recover"));
1500
1501 while (j < 500 && is_ok == 0) {
1502 fa = (*Function)(laststore - j * 1.0e-4);
1503 fb = (*Function)(laststore + j * 1.0e-4);
1504 if ( (fa <= (-FLT_EPSILON) && fb >= FLT_EPSILON)
1505 || (fb <= (-FLT_EPSILON) && fa >= FLT_EPSILON)) {
1506 is_ok = 1;
1507 a = laststore - j * 1e-4;
1508 b = laststore + j * 1e-4;
1509 } else {
1510 ++j;
1511 }
1512 }
1513
1514 /* no success */
1515
1516 if (is_ok == 0) {
1517 if (Flags.interactive == OFF) {
1518 fprintf(stderr,
1519 _("** ERROR **: %s: RootFind: Root not bracketed \n"), caller);
1520 fprintf(stderr,
1521 _("This is likely due to an extreme combination of\n"));
1522 fprintf(stderr, _("mass ratio and fill factor\n"));
1523 abort();
1524 } else {
1525 fprintf(stderr,
1526 _("** ERROR **: %s: RootFind: Root not bracketed \n"), caller);
1527 fprintf(stderr,
1528 _("This is likely due to an extreme combination of\n"));
1529 fprintf(stderr, _("mass ratio and fill factor\n"));
1530 *Err = 1;
1531 return(0.0);
1532 }
1533 }
1534 }
1535
1536 /* -------------- initialize c -------------------------------------- */
1537
1538 fc = fb;
1539
1540 /* -------------- loop until convergence ------------------------------ */
1541
1542 for (i = 0; i < MAXITER; ++i) {
1543
1544 /* If (Root between a and b) then put it between b and c */
1545
1546 if ((fc >= FLT_EPSILON && fb >= FLT_EPSILON)
1547 || (fc <= (-FLT_EPSILON) && fb <= (-FLT_EPSILON) )) {
1548 c = a;
1549 fc = fa;
1550 e = Delta = b-a;
1551 }
1552
1553 /* make fb nearest to x-axis */
1554
1555 if (fabs(fc) <= fabs(fb)) {
1556 a = b; fa = fb;
1557 b = c; fb = fc;
1558 c = a; fc = fa;
1559 }
1560
1561 /* this is the bisection step */
1562
1563 Bisect = (c-b) / 2.0;
1564
1565 /* set tolerance for convergence testing */
1566 /* in the following sum, the first term should ensure proper */
1567 /* convergence testing in the case of a very steep function */
1568
1569 Tolerance1 = (2.0 * FLT_EPSILON * fabs(b)) + (0.5 * Tolerance);
1570
1571 /* Convergence Test */
1572
1573 if ((fabs(Bisect) <= Tolerance1) || (fabs(fb) <= (2.0*FLT_EPSILON))) {
1574 laststore = b;
1575 return(b);
1576 }
1577
1578 /* Interpolation Step if Fast Enough */
1579
1580 if ((fabs(e) >= Tolerance1) && ((fabs(a) - fabs(b)) >= FLT_EPSILON )) {
1581
1582 s = fb/fa;
1583 if (fabs(a - c) <= FLT_EPSILON) {
1584 /* - only two function values, do linear interpolation ------ */
1585 p = 2.0 * Bisect * s;
1586 q = 1.0 - s;
1587 } else {
1588 /* - inverse quadratic interpolation ------- */
1589 q = fa/fc;
1590 r = fb/fc;
1591 p = s * (2.0*Bisect*q*(q-r)-(b-a)*(r-1.0));
1592 q = (q - 1.0) * ( r - 1.0) * (s - 1.0);
1593 }
1594
1595 if (p <= 0.0) p = -p;
1596 else q = -q;
1597
1598 s = e;
1599 e = Delta;
1600
1601 /* Check Bounds */
1602
1603 test1 = (3.0 * Bisect * q) - fabs(Tolerance1 * q);
1604 test2 = fabs( s * q );
1605 if ( ((2.0*p) >= test1) || (p >= test2) ) {
1606 /* we resort to bisection */
1607 Delta = Bisect;
1608 e = Delta;
1609 } else {
1610 /* we accept interpolation */
1611 Delta = p/q;
1612 }
1613
1614 /* Interpolation is creepingly slow, do Bisection */
1615 } else {
1616 Delta = Bisect;
1617 e = Delta;
1618 }
1619
1620 /* we reset a to the former value of b */
1621
1622 a = b; fa = fb;
1623
1624 /* and set b to its updated value */
1625 /* we step at least by Tolerance1 */
1626
1627 if (fabs(Delta) >= Tolerance1) {
1628 b = b + Delta;
1629 } else {
1630 if (Bisect >= DBL_EPSILON) {
1631 Eval = fabs(Tolerance1);
1632 } else {
1633 Eval = -fabs(Tolerance1);
1634 }
1635 b = b + Eval;
1636 }
1637 fb = ((*Function)(b));
1638
1639 }
1640 return(0.0);
1641 }
1642
1643
1644
1645