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 <sys/types.h>
21 #include <math.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include <string.h>
25 #include <float.h>
26
27 #include "Light.h"
28 #include "LightPassband.h"
29
30 #ifdef TM_IN_SYS_TIME
31 #include <sys/time.h>
32 #else
33 #include <time.h>
34 #endif
35
36 #include <pwd.h>
37
38 #ifdef HAVE_UNISTD_H
39 #include <unistd.h>
40 #endif
41
42 #ifdef _WITH_PGPLOT
43 #ifndef _WITH_GNUPLOT
44 #include "cpgplot.h"
45 #endif
46 #endif
47
48
49 /****************************************************************************
50 @package nightfall
51 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
52 @version 1.0
53 @short Initialize Output flux tables
54 @param (void)
55 @return (void)
56 @heading Output
57 ****************************************************************************/
InitOutFlux()58 void InitOutFlux()
59 {
60 int j, band; /* Loop count */
61
62 for (j = 0; j < PhaseSteps; ++j) { /* loop over phase */
63 for (band = 0; band < NUM_MAG; ++band) {
64 FluxOut[j].Mag[band] = 0;
65 }
66 }
67 return;
68 }
69
70
71 /****************************************************************************
72 @package nightfall
73 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
74 @version 1.0
75 @short Write Output flux tables
76 @param (void)
77 @return (void)
78 @heading Output
79 ****************************************************************************/
WriteOutput()80 void WriteOutput()
81 {
82 FILE *file_out; /* output filehandle */
83 unsigned long j; /* Loop count */
84 int band; /* Loop count */
85 double RVp = 0, RVs = 0; /* radial velocities */
86
87 file_out = fopen(OUT_FILE, "w");
88
89 if (file_out == NULL)
90 nf_error(_(errmsg[3]));
91
92 /* ----------- print header ------------------------------------ */
93
94 OutfileHeader(file_out);
95
96
97
98 fprintf(file_out, "# \n");
99 fprintf(file_out, "# -------------------------------------------------\n");
100 fprintf(file_out, "# \n");
101 fprintf(file_out, _("# <Normalization> at <Phase> was:\n"));
102 for (band = 0; band < NUM_MAG; ++band) {
103 fprintf(file_out, _("# Band %3d: %15.6g %15.6g\n"),
104 band, F90[band], P90[band] );
105 }
106
107 fprintf(file_out, "# \n");
108 fprintf(file_out, "# -------------------------------------------------\n");
109 fprintf(file_out, "# \n");
110 fprintf(file_out, _("# SEQ PHASE "));
111 if (Flags.monochrome == OFF)
112 for (band = 0; band < NUM_MAG; ++band) fprintf(file_out, "%16s ", nf_pb_get_name(band));
113 else
114 for (band = 0; band < NUM_MAG; ++band) fprintf(file_out, "%16.4f ", nf_pb_get_wave(band));
115 fprintf(file_out, _("RV(Primary) RV(Secondary) \n"));
116 fprintf(file_out, "# \n");
117
118
119 /* ----------- print data ------------------------------------ */
120
121
122 for (j = 0; j < N_FluxOutFull; ++j) { /* loop over phase */
123
124 if (Flags.elliptic == OFF){
125 RVp = -(FluxOutFull[j].RV[Primary]) /1000.;
126 RVs = -(FluxOutFull[j].RV[Secondary])/1000.;
127 } else {
128 RVp = (FluxOutFull[j].RV[Primary]) /1000.;
129 RVs = (FluxOutFull[j].RV[Secondary])/1000.;
130 }
131
132 /* 2002-05-03 rwichmann: more significant digits
133 */
134 fprintf(file_out, "%6ld %7.4f %16.6f %16.6f %16.6f %16.6f %16.6f %16.6f %16.6f %16.6f %16.6f %16.6f %16.6f %16.6f %11.4g %11.4g\n",
135 j, ((FluxOutFull[j].Phase/(2.0*PI)) - 0.5),
136 FluxOutFull[j].Mag[Umag],
137 FluxOutFull[j].Mag[Bmag],
138 FluxOutFull[j].Mag[Vmag],
139 FluxOutFull[j].Mag[Rmag],
140 FluxOutFull[j].Mag[Imag],
141 FluxOutFull[j].Mag[Jmag],
142 FluxOutFull[j].Mag[Hmag],
143 FluxOutFull[j].Mag[Kmag],
144 FluxOutFull[j].Mag[umag],
145 FluxOutFull[j].Mag[vmag],
146 FluxOutFull[j].Mag[bmag],
147 FluxOutFull[j].Mag[ymag],
148 RVp,
149 RVs);
150
151 }
152
153 (void) fclose(file_out);
154
155 return;
156 }
157
158
159 /****************************************************************************
160 @package nightfall
161 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
162 @version 1.0
163 @short Plot the light curve
164 @param (void)
165 @return (void)
166 @heading Plotting
167 ****************************************************************************/
PlotOutput()168 void PlotOutput()
169 {
170
171 /* >>>>>>>>>>>>>>>>>>>> PLOT <<<<<<<<<<<<<<<<<<<<<<<<< */
172
173
174 #ifdef _WITH_PGPLOT
175
176
177 float MinMag, MaxMag; /* plot scaling */
178 float MinMagP, MaxMagP; /* plot scaling */
179 float MinMagS, MaxMagS; /* plot scaling */
180 float *Xplot, *Yplot, *Zplot; /* data */
181 float *XPplot, *YPplot; /* residuals */
182 int j, j2; /* Loop count */
183 int jstat; /* temp variable */
184 int FlagBand; /* passband to plot */
185 char Title[64]; /* title string */
186 float X1=0.0, X2=1.0; /* ranges */
187 float Y1=-1.0, Y2=0.0; /* ranges */
188 float CalcMax = 0.75; /* upper limit (phase) */
189 char file_out[256+4]; /* ps output file */
190 float xtick, ytick; /* tickmarks */
191 float margin = 0.0; /* above/below data */
192
193 float PLimLow, PLimHigh; /* Limits for minimum */
194
195 int nFactor = 2 * NumPhases;
196 int nOff = NumPhases * PhaseSteps;
197
198 sprintf(file_out, "%s%s", Out_Plot_File, "/CPS");
199
200 PLimHigh = 0.55;
201 PLimLow = 0.45;
202
203 if (Flags.elliptic != OFF) {
204 if ((Orbit.Contact[4] - Orbit.Contact[3]) > DBL_EPSILON) {
205 double cdiff = (Orbit.Contact[4] - Orbit.Contact[3])/3.0;
206 PLimHigh = ((Orbit.Contact[4]+cdiff)/(2.0*PI)) - 0.5;
207 if (PLimHigh < 0.0) PLimHigh += 1.0;
208 PLimLow = ((Orbit.Contact[3]-cdiff)/(2.0*PI)) - 0.5;
209 if (PLimLow < 0.0) PLimLow += 1.0;
210 }
211 }
212
213 /* --------- allocate memory ------------------------- */
214
215 Xplot = malloc(nFactor * PHASESTEPS * sizeof(float));
216 Yplot = malloc(nFactor * PHASESTEPS * sizeof(float));
217 Zplot = malloc(nFactor * PHASESTEPS * sizeof(float));
218
219 XPplot = malloc(4 * MAXOBS * sizeof(float));
220 YPplot = malloc(4 * MAXOBS * sizeof(float));
221
222 /* check allocation */
223
224 if (Xplot == NULL || Yplot == NULL || Zplot == NULL ||
225 XPplot == NULL || YPplot == NULL)
226 #ifdef _WITH_GTK
227 {
228 if ( Flags.interactive == ON) {
229 if (Xplot != NULL) free(Xplot);
230 if (Yplot != NULL) free(Yplot);
231 if (Zplot != NULL) free(Zplot);
232 if (XPplot != NULL) free(XPplot);
233 if (YPplot != NULL) free(YPplot);
234 make_dialog(_(errmsg[0]));
235 return;
236 } else nf_error(_(errmsg[0]));
237 }
238 #else
239 nf_error(_(errmsg[0]));
240 #endif
241
242
243 /* if we plot a lightcurve */
244
245 if (Flags.PlotBand < NUM_MAG) {
246
247 MinMag = 100000.; MaxMag = -100000.;
248 MinMagP = 100000.; MaxMagP = -100000.;
249 MinMagS = 100000.; MaxMagS = -100000.;
250
251 FlagBand = Flags.PlotBand;
252
253 /* ---------- loop over phase ----------------------- */
254
255 for (j = 0; j < nOff/* (NumPhases * PhaseSteps) */; ++j) {
256
257 Xplot[j] = (FluxOutFull[j].Phase/(2.0*PI)) - 0.5;
258 Xplot[j+nOff] = Xplot[j] + NumPhases;
259
260 Yplot[j] = - FluxOutFull[j].Mag[FlagBand];
261 Yplot[j+nOff] = - FluxOutFull[j].Mag[FlagBand];
262
263 MinMag = MIN(MinMag,Yplot[j]); MaxMag = MAX(MaxMag,Yplot[j]);
264 if ((Xplot[j] > PLimLow) && (Xplot[j] < PLimHigh)) {
265 MinMagP = MIN(MinMagP,Yplot[j]); MaxMagP = MAX(MaxMagP,Yplot[j]);
266 }
267 if ((Xplot[j] > -0.1) && (Xplot[j] < 0.1)) {
268 MinMagS = MIN(MinMagS,Yplot[j]); MaxMagS = MAX(MaxMagS,Yplot[j]);
269 }
270
271 }
272
273 CalcMax = Xplot[nOff-1];
274
275 if (Flags.eps == OFF) {
276 if(cpgopen("/XSERVE") != 1)
277 #ifdef _WITH_GTK
278 {
279 if ( Flags.interactive == ON) {
280 free(Xplot);
281 free(Yplot);
282 free(Zplot);
283 free(XPplot);
284 free(YPplot);
285 make_dialog (_(errmsg[2]));
286 return;
287 } else nf_error(_(errmsg[2]));
288 }
289 #else
290 nf_error(_(errmsg[2]));
291 #endif
292 } else {
293 if(cpgopen(file_out) != 1)
294 #ifdef _WITH_GTK
295 {
296 if ( Flags.interactive == ON) {
297 free(Xplot);
298 free(Yplot);
299 free(Zplot);
300 free(XPplot);
301 free(YPplot);
302 make_dialog (_(errmsg[1]));
303 return;
304 } else nf_error(_(errmsg[1]));
305 }
306 #else
307 nf_error(_(errmsg[1]));
308 #endif
309 ++Flags.eps;
310 }
311
312 /* ------------------ start plot ----------------------- */
313
314 #ifdef _WITH_GNUPLOT
315 gnu_start();
316 #endif
317
318 cpgscf(2); cpgslw(1); cpgsch(0.9);
319
320 xtick = 0.0;
321
322 if (Flags.plot == ON) {
323 margin = 0.1 * fabs(MaxMag - MinMag);
324 X1 = -0.25;
325 X2 = (NumPhases - 1.0) + 0.75; /* X2=0.75; */
326 Y1 = MinMag - margin;
327 Y2 = MaxMag + margin;
328 xtick = 0.5;
329 }
330 else if (Flags.plot == 2 ) {
331 margin = 0.1 * fabs(MaxMag - MinMag);
332 X1 = -0.25;
333 X2 = ((2*NumPhases) - 1.0) + 0.75; /* X2=1.75; */
334 Y1 = MinMag - margin;
335 Y2 = MaxMag + margin;
336 xtick = 0.5;
337 }
338 if (Flags.plot == 3 ) {
339 margin = 0.1 * fabs(MaxMag - MinMag);
340 X1= PLimLow; X2=PLimHigh;
341 Y1= MinMagP - margin;
342 Y2= MaxMagP + margin;
343 xtick = (X2-X1)/5.;
344 }
345 if (Flags.plot == 4 ) {
346 margin = 0.1 * fabs(MaxMag - MinMag);
347 X1= -0.05; X2=0.05;
348 Y1= MinMagS - margin;
349 Y2= MaxMagS + margin;
350 xtick = 0.025;
351 }
352
353 ytick = 0.01 * (int)(100 * fabs( (Y1 - Y2) / 6.0 ));
354
355 if (Flags.Passbands[FlagBand] > 0) {
356 cpgsvp(0.2, 0.9, 0.4, 0.9);
357 cpgswin(X1, X2, Y1, Y2);
358 #ifdef _WITH_GNUPLOT
359 cpgbox("BCNST", xtick, 0, "BCNST", ytick, 0);
360 #else
361 cpgbox("BCST", 0, 0, "BCNST", 0, 0);
362 (void) xtick; (void) ytick;
363 #endif
364 } else {
365 cpgsvp(0.2, 0.9, 0.2, 0.9);
366 cpgswin(X1, X2, Y1, Y2);
367 #ifdef _WITH_GNUPLOT
368 cpgbox("BCNST", xtick, 0, "BCNST", ytick, 0);
369 #else
370 cpgbox("BCNST", 0, 0, "BCNST", 0, 0);
371 (void) xtick; (void) ytick;
372 #endif
373 }
374
375 #ifdef _WITH_GNUPLOT
376 cpglab(_("Phase"), _("Brightness (Delta mag)"), "");
377 sprintf(Title, "%s", Filters[FlagBand]);
378 cpgtext(0.0, MaxMag + 0.5*margin, Title);
379 #endif
380
381 /* ---------------------- no data ------------------- */
382
383 if (Flags.Passbands[FlagBand] == 0) {
384
385 cpgline(2*NumPhases*PhaseSteps, Xplot, Yplot);
386 }
387
388 /* ---------------------- observed data ------------- */
389
390 else if (Flags.Passbands[FlagBand] > 0) {
391
392 Y1 = Observed[FlagBand][0].residual;
393 Y2 = Observed[FlagBand][0].residual;
394
395 for (j = 0; j < Flags.Passbands[FlagBand]; ++j)
396 {
397 XPplot[j] = Observed[FlagBand][j].phi;
398 if (XPplot[j] > CalcMax && NumPhases == 1)
399 XPplot[j] = XPplot[j] - 1.0;
400 YPplot[j] = -Observed[FlagBand][j].mag;
401 XPplot[j+Flags.Passbands[FlagBand]] = NumPhases + XPplot[j];
402 YPplot[j+Flags.Passbands[FlagBand]] = YPplot[j];
403 Y1 = MIN(Y1,Observed[FlagBand][j].residual);
404 Y2 = MAX(Y2,Observed[FlagBand][j].residual);
405 }
406
407
408 #ifndef _WITH_GNUPLOT
409 sprintf(Title, "%s", Filters[FlagBand]);
410 cpgtext(-0.1, MaxMag + 0.5*margin, Title);
411
412 cpgline(2 * PhaseSteps * NumPhases, Xplot, Yplot);
413 cpgpt(2*j, XPplot, YPplot, 17);
414 cpglab("", _("Brightness (Delta mag)"), "");
415 #else
416 cpglinept(2 * PhaseSteps * NumPhases, 2*j, Xplot, Yplot, XPplot, YPplot);
417 #endif
418
419 /* ---------------------- residuals ------------------- */
420
421
422 #ifndef _WITH_GNUPLOT
423 cpgsvp(0.2, 0.9, 0.2, 0.39);
424 #else
425 /* no second plot if hardcopy - multiplot not supp. */
426 if (Flags.eps < ON){
427 cpgsvp(0.2, 0.9, 0.0, 0.4);
428 #endif
429 if (Y1 == Y2)
430 {
431 Y1 = -1.0;
432 Y2 = 1.0;
433 }
434 cpgswin(X1, X2, Y1, Y2);
435 ytick = 0.0001 * (int)(10000 * fabs( (Y1 - Y2) / 6));
436 #ifndef _WITH_GNUPLOT
437 cpgbox("ABCNST", 0, 0, "BCNSTV", ytick, 0);
438 /* cpglab("", _("Delta mag"), ""); */
439 #else
440 cpgbox("ABCNST", 0, 0, "BCNST", ytick, 0);
441 cpglab("", _("Delta mag"), "");
442 #endif
443
444 j2 = Flags.Passbands[FlagBand];
445
446 for (j = 0; j < Flags.Passbands[FlagBand]; ++j) {
447 YPplot[j] = Observed[FlagBand][j].residual;
448 YPplot[j2] = Observed[FlagBand][j].residual;
449 ++j2;
450 }
451 cpgpt(2*j, XPplot, YPplot, 17);
452 #ifdef _WITH_GNUPLOT
453 }
454 #endif
455 }
456
457 #ifdef _WITH_GNUPLOT
458 cpglab(_("Phase"), _("Brightness (Delta mag)"), "");
459 #else
460 cpglab(_("Phase"), "", "");
461 #endif
462 if (Flags.eps != OFF) my_cpgend();
463 else cpgend();
464 }
465
466
467 /* >>>>>>>>>>>>>> RADIAL VELOCITY <<<<<<<<<<<<<<<<<<< */
468
469
470 if (Flags.PlotBand == NUM_MAG || Flags.PlotBand == (NUM_MAG+1) ) {
471
472 MinMag = FluxOut[0].RV[Primary]/1000.;
473 MaxMag = FluxOut[0].RV[Primary]/1000.;
474 MinMagP = FluxOut[0].RV[Primary]/1000.;
475 MaxMagP = FluxOut[0].RV[Primary]/1000.;
476 MinMagS = FluxOut[0].RV[Secondary]/1000.;
477 MaxMagS = FluxOut[0].RV[Secondary]/1000.;
478
479 FlagBand = Flags.PlotBand;
480
481 /* -------------- loop over phase -------------------- */
482
483 for (j = 0; j < nOff; ++j) { /* loop over phase */
484
485 Xplot[j] = (FluxOutFull[j].Phase/(PI+PI)) - 0.5;
486 Xplot[j+PhaseSteps] = Xplot[j] + NumPhases;
487
488 if (Flags.elliptic == ON) {
489 Yplot[j] = FluxOutFull[j].RV[Primary]/1000.;
490 Yplot[j+nOff] = FluxOutFull[j].RV[Primary]/1000.;
491 Zplot[j] = FluxOutFull[j].RV[Secondary]/1000.;
492 Zplot[j+nOff] = FluxOutFull[j].RV[Secondary]/1000.;
493 } else {
494 Yplot[j] = - FluxOutFull[j].RV[Primary]/1000.;
495 Yplot[j+nOff] = - FluxOutFull[j].RV[Primary]/1000.;
496 Zplot[j] = - FluxOutFull[j].RV[Secondary]/1000.;
497 Zplot[j+nOff] = - FluxOutFull[j].RV[Secondary]/1000.;
498 }
499
500 MinMagP = MIN(MinMagP,Yplot[j]); MaxMagP = MAX(MaxMagP,Yplot[j]);
501 MinMagS = MIN(MinMagS,Zplot[j]); MaxMagS = MAX(MaxMagS,Zplot[j]);
502
503 }
504
505 CalcMax = Xplot[nOff-1];
506
507 MinMag = MIN(MinMagP,MinMagS);
508 MaxMag = MAX(MaxMagP,MaxMagS);
509
510 /* -------------- open plot -------------------- */
511
512 if (Flags.eps == OFF) {
513 if(cpgopen("/XSERVE") != 1)
514 #ifdef _WITH_GTK
515 {
516 if ( Flags.interactive == ON) {
517 free(Xplot);
518 free(Yplot);
519 free(Zplot);
520 free(XPplot);
521 free(YPplot);
522 make_dialog (_(errmsg[2]));
523 return;
524 } else nf_error(_(errmsg[2]));
525 }
526 #else
527 nf_error(_(errmsg[2]));
528 #endif
529 } else {
530 if(cpgopen(file_out) != 1)
531 #ifdef _WITH_GTK
532 {
533 if ( Flags.interactive == ON) {
534 free(Xplot);
535 free(Yplot);
536 free(Zplot);
537 free(XPplot);
538 free(YPplot);
539 make_dialog (_(errmsg[1]));
540 return;
541 } else nf_error(_(errmsg[1]));
542 }
543 #else
544 nf_error(_(errmsg[1]));
545 #endif
546 ++Flags.eps;
547 }
548
549 /* -------------- plot -------------------- */
550
551 #ifdef _WITH_GNUPLOT
552 gnu_start();
553 #endif
554
555 cpgscf(2); cpgslw(1); cpgsch(0.9);
556
557 xtick = 0.0;
558
559 if (Flags.plot == ON) {
560 margin = 0.1 * fabs(MaxMag - MinMag);
561 X1 = -0.25;
562 X2 = (NumPhases-1) + 0.75;
563 Y1 = MinMag - margin;
564 Y2 = MaxMag + margin;
565 xtick = 0.5;
566 }
567 else if (Flags.plot == 2 ) {
568 margin = 0.1 * fabs(MaxMag - MinMag);
569 X1 = -0.25;
570 X2 = ((2*NumPhases) - 1) + 0.75;
571 Y1 = MinMag - margin;
572 Y2 = MaxMag + margin;
573 xtick = 0.5;
574 }
575 if (Flags.plot == 3 ) {
576 margin = 0.1 * fabs(MaxMag - MinMag);
577 X1= PLimLow; X2=PLimHigh;
578 Y1= MinMagP - margin;
579 Y2= MaxMagP + margin;
580 xtick = (X2-X1)/5.;
581 }
582 if (Flags.plot == 4 ) {
583 margin = 0.1 * fabs(MaxMag - MinMag);
584 X1= -0.1; X2=0.1;
585 Y1= MinMagS - margin;
586 Y2= MaxMagS + margin;
587 xtick = 0.05;
588 }
589
590 if (Flags.Passbands[NUM_MAG] > 0 || Flags.Passbands[NUM_MAG+1] > 0) {
591 if (Y1 == Y2) { Y1 = -20.; Y2 = 20.; } /* a dirty hack */
592 if (Flags.Passbands[NUM_MAG+1] == 0 && Flags.plot == ON) {
593 Y1=MinMagP-0.01; Y2=MaxMagP+0.01;
594 }
595 if (Flags.Passbands[NUM_MAG] == 0 && Flags.plot == ON) {
596 Y1=MinMagS-0.01; Y2=MaxMagS+0.01;
597 }
598 cpgsvp(0.2, 0.9, 0.4, 0.9);
599 cpgswin(X1, X2, Y1, Y2);
600 ytick = 0.1 * (int)(10 * fabs( (Y1 - Y2) / 5));
601 #ifdef _WITH_GNUPLOT
602 cpgbox("BCNST", xtick, 0, "BCNST", ytick, 0);
603 #else
604 cpgbox("BCST", 0, 0, "BCNST", 0, 0);
605 #endif
606 } else {
607 ytick = 0.1 * (int)(10 * fabs( (Y1 - Y2) / 6));
608 cpgsvp(0.2, 0.9, 0.2, 0.9);
609 cpgswin(X1, X2, Y1, Y2);
610 #ifdef _WITH_GNUPLOT
611 cpgbox("BCNST", xtick, 0, "BCNST", ytick, 0);
612 #else
613 cpgbox("BCNST", 0, 0, "BCNST", 0, 0);
614 #endif
615 }
616
617 /* -------------- get data -------------------- */
618
619 if (Flags.Passbands[NUM_MAG] > 0) {
620
621 FlagBand = NUM_MAG;
622
623 Y1 = 1e20;
624 Y2 = -1e20;
625
626 for (j = 0; j < Flags.Passbands[FlagBand]; ++j) {
627
628 XPplot[j] = Observed[FlagBand][j].phi;
629 if (XPplot[j] > CalcMax && NumPhases == 1)
630 XPplot[j] = XPplot[j] - 1.0;
631
632 YPplot[j] = Observed[FlagBand][j].mag ;
633
634 XPplot[j+Flags.Passbands[FlagBand]] = NumPhases + XPplot[j];
635 YPplot[j+Flags.Passbands[FlagBand]] = YPplot[j];
636
637 Y1 = MIN(Y1,Observed[FlagBand][j].residual);
638 Y2 = MAX(Y2,Observed[FlagBand][j].residual);
639
640 }
641 }
642
643 jstat = 2 * Flags.Passbands[NUM_MAG];
644
645 if (Flags.Passbands[NUM_MAG+1] > 0) {
646
647 FlagBand = NUM_MAG+1;
648
649 if (Flags.Passbands[NUM_MAG] > 0) {
650 Y1 = MIN(Y1,Observed[FlagBand][0].residual);
651 Y2 = MAX(Y2,Observed[FlagBand][0].residual);
652 } else {
653 Y1 = 1e20;
654 Y2 = -1e20;
655 }
656
657 /* count up from current state */
658
659 for (j = jstat;
660 j < (jstat + Flags.Passbands[FlagBand]); ++j) {
661
662 XPplot[j] = Observed[FlagBand][j-jstat].phi;
663 if (XPplot[j] > CalcMax && NumPhases == 1)
664 XPplot[j] = XPplot[j] - 1.0;
665
666 YPplot[j] = Observed[FlagBand][j-jstat].mag;
667
668 XPplot[j+Flags.Passbands[FlagBand]] = NumPhases + XPplot[j];
669 YPplot[j+Flags.Passbands[FlagBand]] = YPplot[j];
670
671 Y1 = MIN(Y1,Observed[FlagBand][j-jstat].residual);
672 Y2 = MAX(Y2,Observed[FlagBand][j-jstat].residual);
673
674 }
675
676 }
677
678 /* -------------- plot everything -------------------- */
679
680 #ifdef _WITH_GNUPLOT
681 cpglab(_("Phase"), _("Radial Velocity (km/s)"), "");
682
683 if (Flags.Passbands[NUM_MAG+1] > 0 || Flags.Passbands[NUM_MAG] > 0 ) {
684 cpgline2pt(2*PhaseSteps*NumPhases, 2*PhaseSteps*NumPhases,
685 (2*Flags.Passbands[NUM_MAG] + 2*Flags.Passbands[NUM_MAG+1]),
686 Xplot, Yplot, Xplot, Zplot, XPplot, YPplot);
687 } else {
688 cpgline2(2*PhaseSteps*NumPhases, 2*PhaseSteps*NumPhases,
689 Xplot, Yplot, Xplot, Zplot);
690 }
691 #else
692 cpgline(2*PhaseSteps*NumPhases, Xplot, Yplot);
693 cpgline(2*PhaseSteps*NumPhases, Xplot, Zplot);
694 if (Flags.Passbands[NUM_MAG+1] > 0 || Flags.Passbands[NUM_MAG] > 0 ) {
695 cpgpt( (2*Flags.Passbands[NUM_MAG] + 2*Flags.Passbands[NUM_MAG+1]) ,
696 XPplot, YPplot, 17);
697 }
698 #endif
699
700 cpglab("", _("Radial Velocity (km/s)"), "");
701
702
703 /* ----------- Now plot the residuals -------------- */
704
705
706 if (Flags.Passbands[NUM_MAG] > 0 ) {
707 FlagBand = NUM_MAG;
708 for (j = 0; j < Flags.Passbands[FlagBand]; ++j) {
709 YPplot[j] = Observed[FlagBand][j].residual;
710 YPplot[j+Flags.Passbands[FlagBand]] = Observed[FlagBand][j].residual;
711 }
712 }
713
714
715 if (Flags.Passbands[NUM_MAG+1] > 0 ) {
716 FlagBand = NUM_MAG+1;
717 for (j = jstat; j < (jstat+Flags.Passbands[FlagBand]); ++j) {
718
719 YPplot[j] = Observed[FlagBand][j-jstat].residual;
720 YPplot[j+Flags.Passbands[FlagBand]] =
721 Observed[FlagBand][j-jstat].residual;
722 }
723 }
724
725 #ifndef _WITH_GNUPLOT
726 if (Flags.Passbands[NUM_MAG] > 0 || Flags.Passbands[NUM_MAG+1] > 0) {
727 cpgsvp(0.2, 0.9, 0.2, 0.4);
728 cpgswin(X1, X2, Y1, Y2);
729 cpgbox("ABCNST", 0, 0, "BCNST", 0, 0);
730 cpgpt( (2*Flags.Passbands[NUM_MAG] + 2*Flags.Passbands[NUM_MAG+1]) ,
731 XPplot, YPplot, 17);
732 cpglab(_("Phase"), _("km/s"), "");
733 }
734 #else
735 if (Flags.Passbands[NUM_MAG] > 0 || Flags.Passbands[NUM_MAG+1] > 0) {
736 if (Flags.eps < ON){ /* no second plot if hardcopy */
737 cpgsvp(0.2, 0.9, 0.0, 0.4);
738
739 if (Y1 == Y2)
740 {
741 Y1 = -1.0;
742 Y2 = 1.0;
743 }
744
745 cpgswin(X1, X2, Y1, Y2);
746 ytick = 0.1 * (int)(10 * fabs( (Y1 - Y2) / 6));
747 cpgbox("ABCNST", 0, 0, "BCNST", ytick, 0);
748 cpglab("", _("km/s"), "");
749 cpgpt( (2*Flags.Passbands[NUM_MAG] + 2*Flags.Passbands[NUM_MAG+1]) ,
750 XPplot, YPplot, 17);
751
752 }
753 }
754 #endif
755
756 if (Flags.eps != OFF) my_cpgend();
757 else cpgend();
758 }
759
760
761 free(Xplot);
762 free(Yplot);
763 free(Zplot);
764 free(XPplot);
765 free(YPplot);
766
767 return;
768
769 #endif
770
771 }
772
773 /****************************************************************************
774 @package nightfall
775 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
776 @version 1.37
777 @short Write header of output file
778 @param (FILE*) file_out The output file
779 @return (void)
780 @heading Plotting
781 ****************************************************************************/
OutfileHeader(FILE * file_out)782 void OutfileHeader(FILE *file_out)
783 {
784
785 long j, i; /* loop counters */
786 int band; /* bandpass */
787 int runs = 0, ul = 0, ll = 0; /* runs, upper+lower limit */
788 double Fillout = 0.0; /* Fillout Factor */
789 double ScaleVolume=1.; /* units conversion */
790 double ScaleMass; /* units conversion */
791 double ResMean[NUM_MAG+2]; /* residuals */
792 double ResDev[NUM_MAG+2]; /* residuals */
793 double fratio1 = 1.0; /* async rotation */
794 double fratio2 = 1.0; /* async rotation */
795 double RochePotential; /* roche potential */
796
797 double G_SI = 6.6726e-11; /* gravitational constant */
798
799 char WhoAmI[32] = "Unknown"; /* user name */
800 char *AsciiTime; /* local time */
801 char deftime[] = "66.6"; /* default time */
802
803 /********* UNIX stuff begin *************/
804 /* */
805 /* */
806 time_t time_now;
807 struct tm *time_ptr;
808
809 #ifdef HAVE_UNISTD_H
810 struct passwd *user_passwd;
811 user_passwd = getpwuid(geteuid());
812 strncpy(WhoAmI, user_passwd->pw_gecos, 31);
813 WhoAmI[31] = '\0';
814 #endif
815 (void) (time(&time_now));
816 time_ptr = localtime(&time_now);
817 if (time_ptr != NULL) AsciiTime = asctime(time_ptr);
818 else AsciiTime = &deftime[0];
819
820 /********* UNIX stuff end ***************/
821
822
823 if (Flags.fill == ON)
824 Fillout = (Binary[Primary].RochePot - Binary[Primary].RCLag1)
825 / (Binary[Primary].RCLag2 - Binary[Primary].RCLag1);
826
827 fprintf(file_out, "# -------------------------------------------------\n");
828 fprintf(file_out, "# \n");
829 fprintf(file_out, _("# Thou invoked Nightfall, and it answered thy plea.\n"));
830 fprintf(file_out, "# \n");
831
832 /* for some stupid reason, asctime() inserts a terminating EOL */
833
834 fprintf(file_out, _("# Thy Name is %s, \n"), WhoAmI );
835 fprintf(file_out, _("# and thy time is %s"), AsciiTime);
836 fprintf(file_out, "# \n");
837 fprintf(file_out, "# -------------------------------------------------\n");
838
839
840 fprintf(file_out, "# \n");
841 if (Flags.elliptic == OFF) {
842 fprintf(file_out, _("# System Parameters (circular orbit):\n"));
843 fprintf(file_out, "# \n");
844 } else {
845 fprintf(file_out, _("# System Parameters (elliptic orbit):\n"));
846 fprintf(file_out, "# \n");
847 fprintf(file_out, _("# <Eccentricity> %8.3f (dimensionless)\n"),
848 Orbit.Excentricity);
849 fprintf(file_out, _("# <Omega> %8.3f (degree)\n"),
850 Orbit.Omega);
851 }
852 fprintf(file_out, _("# <Inclination> %8.3f (degree)\n"),
853 (RTOD * Orbit.Inclination));
854 fprintf(file_out, _("# <Mass Ratio> %8.3f (dimensionless)\n"),
855 Binary[Primary].Mq);
856
857 if (Flags.edu != ON) {
858 fprintf(file_out, _("# <Lagrange One> %8.3f (dimensionless)\n"),
859 Binary[Primary].RLag1);
860 }
861
862 /* provide details about the computational model used */
863
864 if (Flags.edu != ON) {
865 fprintf(file_out, "# \n");
866 fprintf(file_out, "# -------------------------------------------------\n");
867 fprintf(file_out, "# \n");
868 fprintf(file_out, _("# Model used:\n"));
869 fprintf(file_out, "# \n");
870 fprintf(file_out, _("# <Flux> %s\n"),
871 Flags.blackbody == ON ? _("Blackbody") : _("Model atmophere"));
872 if (Flags.blackbody == OFF) {
873 fprintf(file_out, _("# <Log g Primary> %3.1f\n"),
874 Binary[Primary].log_g);
875 fprintf(file_out, _("# <Log g Secondary> %3.1f\n"),
876 Binary[Secondary].log_g);
877 }
878 fprintf(file_out, _("# <Fractional visibility> %s\n"),
879 Flags.fractional == ON ? _("On") : _("Off"));
880 fprintf(file_out, _("# <Limb darkening> %s\n"),
881 Flags.limb == 0 ? _("linear") :
882 (Flags.limb == 1 ? _("quadratic") :
883 (Flags.limb == 2 ? _("square root") : _("inverse"))));
884 fprintf(file_out, _("# <Detailed reflection> %d\n"), Flags.reflect);
885 }
886
887 fprintf(file_out, "# \n");
888 fprintf(file_out, "# -------------------------------------------------\n");
889 fprintf(file_out, "# \n");
890 fprintf(file_out, _("# System Size (absolute units):\n"));
891 fprintf(file_out, "# \n");
892
893 if (Flags.edu != ON) {
894
895 fprintf(file_out, _("# <Period> %8.3e (seconds)\n"),
896 Orbit.TruePeriod);
897 fprintf(file_out, _("# <Total Mass> %8.3e (kilogramm)\n"),
898 Orbit.TrueMass);
899 fprintf(file_out, _("# <Distance> %8.3e (meter)\n"),
900 Orbit.TrueDistance);
901 fprintf(file_out, _("# <Lagrange One> %8.3e (meter)\n"),
902 Binary[Primary].RLag1 * PERIDIST);
903 }
904
905 fprintf(file_out, "# \n");
906 fprintf(file_out, _("# <Period> %8.3f (days)\n"),
907 Orbit.TruePeriod / 86400.);
908 fprintf(file_out, _("# <Total Mass> %8.3f (solar mass)\n"),
909 Orbit.TrueMass / 1.989E30);
910 fprintf(file_out, _("# <Distance> %8.3f (solar radius)\n"),
911 Orbit.TrueDistance / SOL_RADIUS);
912 fprintf(file_out, _("# <Lagrange One> %8.3f (solar radius)\n"),
913 Binary[Primary].RLag1 * PERIDIST / SOL_RADIUS);
914
915 fprintf(file_out, "# \n");
916 fprintf(file_out, "# -------------------------------------------------\n");
917 fprintf(file_out, "# \n");
918
919 ScaleMass = Binary[Primary].Mq / ( 1.0 + Binary[Primary].Mq);
920
921 fprintf(file_out, _("# Component Parameters:\n"));
922 fprintf(file_out, _("# Primary Secondary\n"));
923
924 fprintf(file_out, _("# <Max Velocity> %8.3f %8.3f (km/sec)\n"),
925 Binary[Primary].K, Binary[Secondary].K);
926 fprintf(file_out, "# \n");
927
928 ScaleVolume = PERIDIST * PERIDIST * PERIDIST;
929
930 if (Flags.edu != ON) {
931
932 fprintf(file_out, _("# <Mass> %8.3f %8.3f (dimensionless)\n"),
933 1.0 - ScaleMass, ScaleMass);
934 fprintf(file_out, _("# <Gravity> %8.3f %8.3f (dimensionless)\n"),
935 Binary[Primary].Gravity, Binary[Secondary].Gravity);
936
937 /* 2002-05-03 rwichmann: more significant digits
938 */
939 fprintf(file_out, _("# <Polar Radius> %8.3f %8.3f (dimensionless)\n"),
940 (Orbit.Dist*Orbit.MinDist) * Binary[Primary].Radius,
941 (Orbit.Dist*Orbit.MinDist) * Binary[Secondary].Radius);
942
943 if (Flags.fill == OFF)
944 {
945 fprintf(file_out,
946 _("# <Point Radius> %8.6f %8.6f (dimensionless)\n"),
947 (Orbit.Dist*Orbit.MinDist) * Binary[Primary].Xp,
948 (Orbit.Dist*Orbit.MinDist) * Binary[Secondary].Xp);
949 fprintf(file_out,
950 _("# <Mean Radius> %8.6f %8.6f (dimensionless)\n"),
951 Orbit.MinDist * Binary[Primary].RadMean,
952 Orbit.MinDist * Binary[Secondary].RadMean);
953 fprintf(file_out,
954 _("# <Volume> %8.6f %8.6f (dimensionless)\n"),
955 Binary[Primary].Volume * Orbit.MinDist * Orbit.MinDist * Orbit.MinDist,
956 Binary[Secondary].Volume * Orbit.MinDist * Orbit.MinDist * Orbit.MinDist);
957 }
958
959 fprintf(file_out, "# \n");
960
961 fprintf(file_out, _("# <Mass> %8.2e %8.2e (kilogramm)\n"),
962 Orbit.TrueMass * (1.0 - ScaleMass), Orbit.TrueMass * ScaleMass);
963
964 fprintf(file_out, _("# <Gravity> %8.3f %8.3f (m kg/s^2)\n"),
965 Binary[Primary].Gravity
966 * Orbit.TrueMass * (1.0 - ScaleMass)
967 * G_SI
968 / (Orbit.TrueDistance * Orbit.TrueDistance),
969 Binary[Secondary].Gravity
970 * Orbit.TrueMass * (1.0 - ScaleMass)
971 * G_SI
972 / (Orbit.TrueDistance * Orbit.TrueDistance));
973
974 fprintf(file_out, _("# <Polar Radius> %8.2e %8.2e (meter)\n"),
975 Binary[Primary].Radius * PERIDIST * Orbit.Dist,
976 Binary[Secondary].Radius * PERIDIST * Orbit.Dist);
977
978 if (Flags.fill == OFF)
979 {
980 fprintf(file_out, _("# <Point Radius> %8.2e %8.2e (meter)\n"),
981 Binary[Primary].Xp * PERIDIST * Orbit.Dist,
982 Binary[Secondary].Xp * PERIDIST * Orbit.Dist);
983 fprintf(file_out, _("# <Mean Radius> %8.2e %8.2e (meter)\n"),
984 Binary[Primary].RadMean * PERIDIST,
985 Binary[Secondary].RadMean * PERIDIST);
986 fprintf(file_out,
987 _("# <Volume> %8.2e %8.2e (cubic meter)\n"),
988 Binary[Primary].Volume * ScaleVolume,
989 Binary[Secondary].Volume * ScaleVolume);
990 }
991 }
992
993 fprintf(file_out, "# \n");
994 fprintf(file_out, _("# <Mass> %8.3f %8.3f (solar mass)\n"),
995 Orbit.TrueMass * (1.0 - ScaleMass) / 1.989E30,
996 Orbit.TrueMass * ScaleMass / 1.989E30 );
997 fprintf(file_out, _("# <Polar Radius> %8.3f %8.3f (solar radius)\n"),
998 Binary[Primary].Radius * PERIDIST * Orbit.Dist / SOL_RADIUS,
999 Binary[Secondary].Radius * PERIDIST * Orbit.Dist / SOL_RADIUS);
1000 if (Flags.fill == OFF) {
1001 fprintf(file_out, _("# <Point Radius> %8.3f %8.3f (solar radius)\n"),
1002 Binary[Primary].Xp * PERIDIST * Orbit.Dist / SOL_RADIUS,
1003 Binary[Secondary].Xp * PERIDIST * Orbit.Dist / SOL_RADIUS);
1004 fprintf(file_out, _("# <Mean Radius> %8.3f %8.3f (solar radius)\n"),
1005 Binary[Primary].RadMean * PERIDIST / SOL_RADIUS,
1006 Binary[Secondary].RadMean * PERIDIST / SOL_RADIUS);
1007 ScaleVolume = ScaleVolume / (SOL_RADIUS * SOL_RADIUS * SOL_RADIUS);
1008 fprintf(file_out, _("# <Volume> %8.3f %8.3f (solar volumes)\n"),
1009 Binary[Primary].Volume * ScaleVolume,
1010 Binary[Secondary].Volume * ScaleVolume);
1011 }
1012 fprintf(file_out, "# \n");
1013 fprintf(file_out, _("# <Mean Temp.> %8.1f %8.1f (Kelvin)\n"),
1014 Binary[Primary].TempMean, Binary[Secondary].TempMean);
1015
1016
1017 fprintf(file_out, "# \n");
1018 fprintf(file_out, "# -------------------------------------------------\n");
1019 fprintf(file_out, "# \n");
1020
1021
1022 fprintf(file_out, _("# Primary Secondary\n"));
1023 if (Flags.fill == OFF) {
1024 fprintf(file_out, _("# <Fill Factor> %8.3f %8.3f (dimensionless)\n"),
1025 Binary[Primary].RocheFill, Binary[Secondary].RocheFill);
1026 } else {
1027 fprintf(file_out, _("# <Fill Factor> %8.3f (dimensionless)\n"),
1028 Binary[Primary].RocheFill);
1029 }
1030 fprintf(file_out, _("# <Temperature> %8.1f %8.1f (Kelvin)\n"),
1031 Binary[Primary].Temperature,
1032 Binary[Secondary].Temperature);
1033
1034 if (Flags.edu != ON) {
1035
1036 if (Flags.asynchron1 == ON) fratio1 = Binary[Primary].Fratio;
1037 if (Flags.asynchron2 == ON) fratio2 = Binary[Secondary].Fratio;
1038 fprintf(file_out, _("# <Asynchronicity>%8.3f %8.3f (dimensionless)\n"),
1039 fratio1, fratio2);
1040
1041 fprintf(file_out, _("# <Albedo> %8.3f %8.3f (dimensionless)\n"),
1042 Binary[Primary].Albedo, Binary[Secondary].Albedo);
1043 fprintf(file_out, _("# <Grav. Dark.> %8.3f %8.3f (dimensionless)\n"),
1044 Binary[Primary].GravDarkCoeff ,
1045 Binary[Secondary].GravDarkCoeff);
1046
1047 fprintf(file_out, "# \n");
1048 fprintf(file_out, "# -------------------------------------------------\n");
1049 fprintf(file_out, "# \n");
1050
1051 fprintf(file_out,
1052 _("# <Surface Potential Primary> %10.4f (dimensionless)\n"),
1053 Binary[Primary].RochePot);
1054 fprintf(file_out,
1055 _("# <Lagrange 1 Potential> %10.4f (dimensionless)\n"),
1056 Binary[Primary].RCLag1);
1057 fprintf(file_out,
1058 _("# <Lagrange 2 Potential> %10.4f (dimensionless)\n"),
1059 Binary[Primary].RCLag2);
1060 if (Flags.fill == ON) {
1061 fprintf(file_out,
1062 _("# <Fillout Factor> %10.4f (dimensionless)\n"),
1063 Fillout);
1064 }
1065 fprintf(file_out,
1066 _("# <Surface Potential Secondary> %10.4f (dimensionless)\n"),
1067 Binary[Secondary].RochePot);
1068
1069 RochePot = 0.0; /* Binary[Primary].RochePot; */
1070 Mq = Binary[Primary].Mq;
1071
1072 RochePotential = RochePerpendSecond(Binary[Secondary].Radius);
1073 /*
1074 RochePotential = 1.0/Binary[Secondary].Radius
1075 + Binary[Secondary].Mq/sqrt(1 + SQR(Binary[Secondary].Radius));
1076 */
1077
1078 fprintf(file_out,
1079 _("# <... in Primary frame> %10.4f (dimensionless)\n"),
1080 RochePotential);
1081
1082 /* Kontaktzeiten */
1083
1084 if ((Orbit.Contact[2] - Orbit.Contact[1]) > DBL_EPSILON)
1085 {
1086 double tdiff = (Orbit.TruePeriod/60.0)/PhaseSteps;
1087 fprintf(file_out,
1088 _("# <Secondary eclipse begin> %10.4f (minutes +/- %5.2f)\n"),
1089 ((Orbit.TruePeriod/60.0) * Orbit.Contact[1]/(2.0*PI)) - tdiff,
1090 tdiff/2.0);
1091 fprintf(file_out,
1092 _("# <Secondary eclipse end> %10.4f (minutes +/- %5.2f)\n"),
1093 ((Orbit.TruePeriod/60.0) * Orbit.Contact[2]/(2.0*PI)) - tdiff,
1094 tdiff/2.0);
1095 }
1096 else
1097 fprintf(file_out,
1098 _("# <Secondary eclipse begin> --no eclipse--\n"));
1099
1100
1101 if ((Orbit.Contact[6] - Orbit.Contact[5]) > DBL_EPSILON)
1102 {
1103 double tdiff = (Orbit.TruePeriod/60.0)/PhaseSteps;
1104 fprintf(file_out,
1105 _("# <Secondary totality begin> %10.4f (minutes +/- %5.2f)\n"),
1106 ((Orbit.TruePeriod/60.0) * Orbit.Contact[5]/(2.0*PI)) - tdiff,
1107 tdiff/2.0);
1108 fprintf(file_out,
1109 _("# <Secondary totality end> %10.4f (minutes +/- %5.2f)\n"),
1110 ((Orbit.TruePeriod/60.0) * Orbit.Contact[6]/(2.0*PI)) - tdiff,
1111 tdiff/2.0);
1112 }
1113 else
1114 fprintf(file_out,
1115 _("# <Secondary totality begin> --no totality--\n"));
1116
1117 if ((Orbit.Contact[4] - Orbit.Contact[3]) > DBL_EPSILON)
1118 {
1119 double tdiff = (Orbit.TruePeriod/60.0)/PhaseSteps;
1120 fprintf(file_out,
1121 _("# <Primary eclipse begin> %10.4f (minutes +/- %5.2f)\n"),
1122 ((Orbit.TruePeriod/60.0) * Orbit.Contact[3]/(2.0*PI)) - tdiff,
1123 tdiff/2.0);
1124 fprintf(file_out,
1125 _("# <Primary eclipse end> %10.4f (minutes +/- %5.2f)\n"),
1126 ((Orbit.TruePeriod/60.0) * Orbit.Contact[4]/(2.0*PI)) - tdiff,
1127 tdiff/2.0);
1128 }
1129 else
1130 fprintf(file_out,
1131 _("# <Primary eclipse begin> --no eclipse--\n"));
1132
1133
1134 if ((Orbit.Contact[8] - Orbit.Contact[7]) > DBL_EPSILON)
1135 {
1136 double tdiff = (Orbit.TruePeriod/60.0)/PhaseSteps;
1137 fprintf(file_out,
1138 _("# <Primary totality begin> %10.4f (minutes +/- %5.2f)\n"),
1139 ((Orbit.TruePeriod/60.0) * Orbit.Contact[7]/(2.0*PI)) - tdiff,
1140 tdiff/2.0);
1141 fprintf(file_out,
1142 _("# <Primary totality end> %10.4f (minutes +/- %5.2f)\n"),
1143 ((Orbit.TruePeriod/60.0) * Orbit.Contact[8]/(2.0*PI)) - tdiff,
1144 tdiff/2.0);
1145 }
1146 else
1147 fprintf(file_out,
1148 _("# <Primary totality begin> --no totality--\n"));
1149
1150
1151
1152 fprintf(file_out, "# \n");
1153 fprintf(file_out, "# -------------------------------------------------\n");
1154 fprintf(file_out, "# \n");
1155 fprintf(file_out, _("# <Luminosities>:\n"));
1156 for (band = 0; band < NUM_MAG; ++band) {
1157 fprintf(file_out, _("# Band %3d: %15.6g %15.6g"),
1158 band, L90[Primary][band], L90[Secondary][band] );
1159 #ifdef HAVE_DISK
1160 if (Flags.disk == ON)
1161 {
1162 fprintf(file_out, " %15.6g\n", L90[Disk][band]);
1163 }
1164 else
1165 {
1166 fprintf(file_out, "\n");
1167 }
1168 #else
1169 fprintf(file_out, "\n");
1170 #endif
1171 }
1172 }
1173
1174 fprintf(file_out, "# \n");
1175 fprintf(file_out, "# -------------------------------------------------\n");
1176 fprintf(file_out, "# \n");
1177
1178 fprintf(file_out, _("# Spots (if any)\n"));
1179 fprintf(file_out, "# \n");
1180 if (Flags.Spots1 > 0) {
1181 fprintf(file_out, "# \n");
1182 fprintf(file_out, _("# <Spots on Primary>\n"));
1183 fprintf(file_out,
1184 _("# Longitude Latitude Radius DimFactor\n"));
1185 for (j = 0; j < Flags.Spots1; ++j) {
1186 fprintf(file_out, _("# <Spot %3ld> %8.3f %8.3f %8.3f %8.3f\n"),
1187 j,
1188 Spot[Primary][j].longitude,
1189 Spot[Primary][j].latitude,
1190 Spot[Primary][j].radius,
1191 Spot[Primary][j].dimfactor );
1192 }
1193 fprintf(file_out, "# \n");
1194 }
1195 if (Flags.Spots2 > 0) {
1196 fprintf(file_out, "# \n");
1197 fprintf(file_out, _("# <Spots on Secondary>\n"));
1198 fprintf(file_out,
1199 _("# Longitude Latitude Radius DimFactor\n"));
1200 for (j = 0; j < Flags.Spots2; ++j) {
1201 fprintf(file_out, _("# <Spot %3ld> %8.3f %8.3f %8.3f %8.3f\n"),
1202 j,
1203 Spot[Secondary][j].longitude,
1204 Spot[Secondary][j].latitude,
1205 Spot[Secondary][j].radius,
1206 Spot[Secondary][j].dimfactor );
1207 }
1208 }
1209
1210 if (Flags.edu != ON) {
1211 fprintf(file_out, "# \n");
1212 fprintf(file_out, "# -------------------------------------------------\n");
1213 fprintf(file_out, "# \n");
1214
1215 fprintf(file_out, _("# Fit Results (if any)\n"));
1216 fprintf(file_out, "# \n");
1217
1218 if (Flags.IsComputed == ON) {
1219
1220 for (band=0; band < (NUM_MAG+2); ++band) {
1221
1222 if (Flags.Passbands[band] > 0) {
1223
1224 Runs(Observed[band], Flags.Passbands[band], &runs, &ul, &ll);
1225
1226 ResMean[band] = 0.0; ResDev[band] = 0.0;
1227 for (i = 0; i < Flags.Passbands[band]; ++i) {
1228 ResMean[band] =
1229 ResMean[band] + Observed[band][i].residual;
1230 }
1231 ResMean[band] = ResMean[band]/Flags.Passbands[band];
1232 for (i = 0; i < Flags.Passbands[band]; ++i) {
1233 ResDev[band] = ResDev[band] +
1234 SQR(ResMean[band]-Observed[band][i].residual);
1235 }
1236 ResDev[band] = sqrt(ResDev[band])/(Flags.Passbands[band]-1);
1237
1238 fprintf(file_out,
1239 _("# <%12s> Mean Residual %9.5g SDV Residuals %9.5f\n"),
1240 Filters[band], ResMean[band], ResDev[band]);
1241
1242 fprintf(file_out,
1243 _("# Runs %5d Upper Limit %5d Lower Limit %6d\n"),
1244 runs, ul, ll);
1245 fprintf(file_out,"# \n");
1246 }
1247 }
1248 }
1249 }
1250
1251 /*@i@*/ return;
1252
1253 }
1254