/* NIGHTFALL Light Curve Synthesis Program */
/* Copyright (C) 1998 Rainer Wichmann */
/* */
/* This program is free software; you can redistribute it */
/* and/or modify */
/* it under the terms of the GNU General Public License as */
/* published by */
/* the Free Software Foundation; either version 2 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
/* GNU General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program; if not, write to the Free Software */
/* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
/* Define this if you want to print coordinates
*/
/* #define PRINT_COORDS 1 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>
#include "Light.h"
#ifdef _WITH_PGPLOT
#ifndef _WITH_GNUPLOT
#include "cpgplot.h"
#endif
#endif
#ifdef _WITH_PGPLOT
/******************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Real-time animation of lightcurve computation
@param (int) j The phase index
@return (void)
@heading Animation
*******************************************************************/
void Animate(int j)
{
long i; /* loop variable */
int test; /* test variable */
int Comp; /* the stellar component */
int MaxComp = 2; /* the maximum # of components */
long CountElem; /* # of elements in X/Yplot */
float *Xplot; /* plot array */
float *Yplot; /* plot array */
float *Zplot; /* plot array */
double sinAlpha; /* sin inclination */
double cosAlpha; /* cos inclination */
double sinPhi; /* sin phase */
double cosPhi; /* cos phase */
double Trans[4][4]; /* transformation matrix */
double t10, t11, t12, t13; /* transformation matrix */
double t20, t21, t22, t23; /* transformation matrix */
double Calib; /* normalization for flux */
double C; /* sign (for primary/secondary) */
float x1, x2, y1, y2; /* plot window borders */
float MinMag, MaxMag; /* plot scaling */
float MinVel, MaxVel; /* plot scaling */
float Upper; /* plot window borders */
float Xtick; /* tickmark distance */
PhotoOut *FluxPtr = FluxOut; /* pointer to flux */
SurfaceElement *SurfPtr; /* pointer to surface */
OrbitType *OrbPtr = &Orbit; /* pointer to Orbit */
char titlestring[MAX_CFG_INLINE+1]; /* the title string */
float com; /* centre of mass */
int js; /* j minus invalid */
#ifdef HAVE_DISK
static float init_phase; /* initial phase */
#endif
#if defined(PRINT_COORDS)
static FILE * outfile;
double sx[3], sy[3];
double ssin, scos, posw, hyp;
#endif
/* >>>>>>>>>>> Initialize <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
Xplot = malloc(3 * MAXELEMENTS * sizeof(float));
Yplot = malloc(3 * MAXELEMENTS * sizeof(float));
Zplot = malloc(PHASESTEPS * sizeof(float));
if (Xplot == NULL || Yplot == NULL || Zplot == NULL)
#ifdef _WITH_GTK
{
if (Flags.interactive == ON) {
make_dialog (_(errmsg[0]));
if (Xplot != NULL) free(Xplot);
if (Yplot != NULL) free(Yplot);
if (Zplot != NULL) free(Zplot);
return;
} else nf_error(_(errmsg[5]));
}
#else
nf_error(_(errmsg[5]));
#endif
sinAlpha = sin(OrbPtr->Inclination - 0.5*PI);
cosAlpha = cos(OrbPtr->Inclination - 0.5*PI);
/* >>>>>>>>>>> open device <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
#ifdef _WITH_GNUPLOT
if (j == 0) {
if (cpgopen("/XSERVE") != 1)
#ifdef _WITH_GTK
{
if (Flags.interactive == ON) {
make_dialog (_(errmsg[2]));
free(Xplot);
free(Yplot);
free(Zplot);
return;
} else nf_error(_(errmsg[2]));
}
#else
nf_error(_(errmsg[2]));
#endif
}
gnu_start();
#else
if(cpgopen("/XSERVE") != 1)
#ifdef _WITH_GTK
{
if (Flags.interactive == ON) {
make_dialog (_(errmsg[2]));
free(Xplot);
free(Yplot);
free(Zplot);
return;
} else nf_error(_(errmsg[2]));
}
#else
nf_error(_(errmsg[2]));
#endif
#endif
/* plot box */
if (Orbit.Name[0] != '\0')
strncpy(titlestring, Orbit.Name, MAX_CFG_INLINE);
else
strcpy(titlestring, "Nightfall");
cpgscf(2); cpgslw(1); cpgsch(0.9);
cpgsvp(0.05,0.6,0.05,0.9);
cpgwnad(-0.98, 0.980, -0.9, 0.9);
/* cpgbox("BCNST", 0,0, "BCNST", 0,0); */
cpgsch(2.0); cpgmtxt("T", 1, 0.8, 0.5, titlestring); cpgsch(0.9);
/* >>>>>>>>>>>> Light Curve <<<<<<<<<<<<<<<<<<<<<<<<<< */
/* photometric zeropoint */
Calib = FluxOut[0].Mag[Vmag];
/* get photometric data and scale */
js = 0;
MinMag = 10.;
MaxMag = -10.;
FluxPtr = FluxOut;
for (i = 0; i <= j; ++i) {
if (FluxPtr->Mag[Vmag] != 0)
{
Xplot[js] = (FluxPtr->Phase/(PI+PI)) - 0.5;
Yplot[js] = 2.5 * log10(FluxPtr->Mag[Vmag]/Calib);
MinMag = MIN(MinMag,Yplot[js]); MaxMag = MAX(MaxMag,Yplot[js]);
++js;
}
++FluxPtr;
}
/* shift phase if necessary */
test = ON;
do {
if (Xplot[0] > 0.0) {
for (i = 0; i <= js; ++i) {
Xplot[i] = Xplot[i] - 1.0;
}
test = ON; /* there was something to do */
} else {
if (Xplot[0] < (-1.0)) {
for (i = 0; i <= js; ++i) {
Xplot[i] = Xplot[i] + 1.0;
}
test = ON; /* there was something to do */
} else { test = OFF; }
}
} while (test == ON);
if (fabs(MinMag-MaxMag) < 0.1) MinMag = -0.1;
/* get lower window border */
cpgqvp(0, &x1, &x2, &y1, &y2);
/* store upper window border */
Upper = y2;
/* plot box */
Xtick = (int)(100.*(MaxMag-MinMag+0.1)/3.)/100.;
cpgscf(2); cpgslw(1); cpgsch(0.9);
cpgsvp(0.6,0.95,y1,y1+(y2-y1)/2.);
/* was commented out */
cpgswin(-0.4, 0.9, MinMag-0.05, MaxMag+0.05);
cpgbox("BCNST", 0.5, 5, "BCMST", Xtick, 4);
cpgmtxt("B", 3, 0.5, 0.5, "Phase");
/* reset point size and plot lightcurve */
cpgsch(0.1); cpgslw(1);
cpgsci(2);
cpgline(js, Xplot, Yplot);
cpgsci(1);
/* >>>>>>>>>>>> Radial Velocities <<<<<<<<<<<<<<<<<< */
FluxPtr = FluxOut;
MinVel = (FluxPtr->RV[Primary]) /1000.;
MaxVel = (FluxPtr->RV[Secondary])/1000.;
js = 0;
if (Flags.elliptic == ON) {
for (i = 0; i <= j; ++i) {
if (FluxPtr->Mag[Vmag] != 0)
{
Yplot[js] = (FluxPtr->RV[Primary]) /1000.;
Zplot[js] = (FluxPtr->RV[Secondary])/1000.;
MinVel = MIN(MinVel,Yplot[js]); MaxVel = MAX(MaxVel,Yplot[js]);
MinVel = MIN(MinVel,Zplot[js]); MaxVel = MAX(MaxVel,Zplot[js]);
++js;
}
++FluxPtr;
}
} else {
for (i = 0; i <= j; ++i) {
if (FluxPtr->Mag[Vmag] != 0)
{
Yplot[js] = - (FluxPtr->RV[Primary]) /1000.;
Zplot[js] = - (FluxPtr->RV[Secondary])/1000.;
MinVel = MIN(MinVel,Yplot[js]); MaxVel = MAX(MaxVel,Yplot[js]);
MinVel = MIN(MinVel,Zplot[js]); MaxVel = MAX(MaxVel,Zplot[js]);
++js;
}
++FluxPtr;
}
}
if (fabs(MinVel-MaxVel) < 5.) {
MinVel = -5.01;
MaxVel = 5.01;
}
/* get lower window border */
cpgqvp(0, &x1, &x2, &y1, &y2);
/* plot box */
cpgscf(2); cpgslw(1); cpgsch(0.9);
cpgsvp(0.6,0.95,y2,Upper);
cpgswin(-0.4, 0.9, MinVel-1., MaxVel+1.);
#ifndef _WITH_GNUPLOT
cpgbox("BCST", 0.5, 5, "BCMST", (int)((2.+MaxVel-MinVel)/3), 4);
#else
cpgbox("BCNST", 0.5, 5, "BCMST", (int)((2.+MaxVel-MinVel)/3), 4);
#endif
/* reset point size and plot lightcurve */
cpgsch(0.1); cpgslw(1);
#ifndef _WITH_GNUPLOT
/* primary */
cpgsci(5); cpgline(js, Xplot, Yplot); cpgsci(1);
/* secondary */
cpgsci(4); cpgline(js, Xplot, Zplot); cpgsci(1);
#else
cpgline2(js, js, Xplot, Yplot, Xplot, Zplot);
#endif
#ifndef _WITH_GNUPLOT
/* plot the non-keplerian disturbance */
/* zoomed by factor two */
js = 0;
FluxPtr = FluxOut;
if (Flags.elliptic == OFF){
for (i = 0; i <= j; ++i) {
if (FluxPtr->Mag[Vmag] != 0)
{
Yplot[js] = -(FluxPtr->D_RV[Primary]) /500.;
Zplot[js] = -(FluxPtr->D_RV[Secondary])/500.;
++js;
}
++FluxPtr;
}
} else {
for (i = 0; i <= j; ++i) {
if (FluxPtr->Mag[Vmag] != 0)
{
Yplot[js] = (FluxPtr->D_RV[Primary]) /500.;
Zplot[js] = (FluxPtr->D_RV[Secondary])/500.;
++js;
}
++FluxPtr;
}
}
cpgsls(4);
/* primary */
cpgsci(5); cpgline(js, Xplot, Yplot); cpgsci(1);
/* secondary */
cpgsci(4); cpgline(js, Xplot, Zplot); cpgsci(1);
cpgsls(1);
#endif
/* >>>>>>>>>>>>> the stars <<<<<<<<<<<<<<<<<<<<<<< */
/* plot box */
cpgscf(2); cpgslw(1); cpgsch(0.9);
cpgsvp(0.05,0.6,0.05,0.9);
cpgwnad(-0.98, 0.980, -0.9, 0.9);
cpgsch(2.0); cpgmtxt("T", 1, 0.8, 0.5, titlestring); cpgsch(0.9);
/* reset size */
cpgsch(0.1); cpgslw(1);
#ifdef HAVE_DISK
if (j == 0) {
/* die scheibe muss raumfest sein, z.B. wg. tilt */
init_phase = OrbPtr->Phase;
}
#endif
/* loop over components */
#if defined(PRINT_COORDS)
if (j == 0)
{
outfile = fopen ("NightfallCoordinates.txt", "w");
if (!outfile)
perror("Error opening NightfallCoordinates.txt");
}
#endif
CountElem = 0;
#ifdef HAVE_DISK
if (Flags.disk == ON)
MaxComp = 3;
#endif
for (Comp = 0; Comp < MaxComp; ++Comp) {
/* now we need the phase */
/* we get it from Orbit.Phase, which is set either */
/* in KeplerSolve or in LightCurve */
/* so this routine must be called after LightCurve */
if (Comp == Primary) {
sinPhi = sin(-OrbPtr->Phase);
cosPhi = cos(-OrbPtr->Phase);
C = 1.0;
} else if (Comp == Secondary) {
C = -1.0; /* mirror y -> (-y) for secondary */
sinPhi = sin(-OrbPtr->Phase + PI);
cosPhi = cos(-OrbPtr->Phase + PI);
} else { /* (Comp == Disk) */
C = -1.0; /* mirror y -> (-y) for secondary */
sinPhi = sin(-OrbPtr->Phase + PI);
cosPhi = cos(-OrbPtr->Phase + PI);
}
/* the centre of mass */
if (Comp == Primary || Comp == Secondary)
com = (((float)(Binary[Comp].Mq))/( 1.0 + (float)(Binary[Comp].Mq)));
else
com = (((float)(Binary[Secondary].Mq))/( 1.0 + (float)(Binary[Secondary].Mq)));
/* the transformation matrix */
/* first column, then row */
/* Trans[0][0] = cosPhi * cosAlpha; */
Trans[1][0] = sinPhi;
Trans[2][0] = cosPhi * sinAlpha;
/* Trans[3][0] = 0.0; */
/* Trans[0][1] = -sinPhi * cosAlpha;*/
Trans[1][1] = cosPhi;
Trans[2][1] = -sinPhi * sinAlpha;
/* Trans[3][1] = 0.0; */
/* Trans[0][2] = -sinAlpha; */
Trans[1][2] = 0.0;
Trans[2][2] = cosAlpha;
/* Trans[3][2] = 0.0; */
/* Trans[0][3] = com * (1.0 - cosPhi * cosAlpha); */
Trans[1][3] = -com * sinPhi;
Trans[2][3] = -com * cosPhi * sinAlpha;
/* Trans[3][3] = 1.0; */
/*
Trans[0][3] = Binary[Comp].RLag1 * (1.0 - cosPhi * cosAlpha);
Trans[1][3] = -Binary[Comp].RLag1 * sinPhi;
Trans[2][3] = -Binary[Comp].RLag1 * cosPhi * sinAlpha;
Trans[3][3] = 1.0;
*/
t10 = Trans[1][0]; t11 = Trans[1][1]; t12 = Trans[1][2]; t13 = Trans[1][3];
t20 = Trans[2][0]; t21 = Trans[2][1]; t22 = Trans[2][2]; t23 = Trans[2][3];
/* do the transformation */
SurfPtr = Surface[Comp];
for (i = 0; i < Binary[Comp].NumElem; ++i) {
if (SurfPtr->EclipseFlag == 0 &&
SurfPtr->SpotNum == 0)
{
/* Observer looks from positive X-Axis */
/* this is the negative Y-Axis */
Xplot[CountElem] = - ( SurfPtr->lambda * t10
+ SurfPtr->mu * t11 * C
+ SurfPtr->nu * t12
+ t13);
/* + 1.0 * t13); */
/* and this the Z-Axis */
Yplot[CountElem] = ( SurfPtr->lambda * t20
+ SurfPtr->mu * t21 * C
+ SurfPtr->nu * t22
+ t23);
/* + 1.0 * t23); */
/* increase counter for elements in Xplot/Yplot */
++CountElem;
}
++SurfPtr;
}
}
#if defined(PRINT_COORDS)
hyp = sqrt(((sx[0] - sx[1])*(sx[0] - sx[1]))+((sy[0] - sy[1])*(sy[0] - sy[1])));
ssin = (sy[0] - sy[1])/hyp;
scos = (sx[0] - sx[1])/hyp;
posw = atan((sy[0] - sy[1])/(sx[0] - sx[1]))*(180/M_PI);
if (ssin < 0 && scos > 0) posw += 180;
if (ssin > 0 && scos > 0) posw += 180;
if (ssin > 0 && scos < 0) posw += 360;
if (outfile)
{
fprintf (outfile, "%3d %12.6f %12.6f %12.6f\n",
j, Orbit.Dist*hyp, posw /* , ssin, scos */,
(FluxOut[j].Phase/(PI+PI)) - 0.5);
fflush(outfile);
}
#endif
/* plot the data */
#ifdef _WITH_GNUPLOT
cpgsvp(0.05,0.6,0.05,0.9);
cpgwnad(-0.98, 0.980, -0.9, 0.9);
cpgbox("BCST", 0, 0, "BCST", 0, 0);
cpgmtxt("T", 1, 0.8, 0.5, titlestring);
#endif
cpgpt(CountElem, Xplot, Yplot, 17);
/* plot finished */
#ifndef _WITH_GNUPLOT
cpgend();
#else
gnu_end();
if (j == (PhaseSteps - 1)) cpgend();
#endif
free(Xplot);
free(Yplot);
free(Zplot);
return;
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1