/* 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. */
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <float.h>
#include "Light.h"
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short The radial velocity curve and the line profile
@param (int) Phase The phase index
@return (void)
@heading Light Curve
****************************************************************************/
void RadVel(int Phase)
{
long i; /* loop counter */
int Comp; /* stellar component */
double Profile[2][PROFILE_ARRAY]; /* profile value */
double Lambda[PROFILE_ARRAY]; /* wavelength */
double Offset, ConstantFactor; /* conversion wl -> index */
int HashLo, HashHi; /* hash variable */
double SpeedOfLight = 2.9979e8; /* in meter/second */
double Cosgamma, Common; /* flux computing */
double OneMinCosgamma, SumOfWeights;/* flux computing */
double Velocity; /* the velocity */
double Weight = 0.0; /* weight */
double WeightLo, WeightHi, Lo, Hi; /* for weight computation */
FILE *file_out; /* output filehandle */
char Out_File[32]; /* output filename */
double CorrPhase; /* phase */
double DeltaLambda; /* step in wavelength */
double axis_1; /* semimajor axis */
double temp; /* temporary */
int VelMag = Vmag; /* bandpass for velocity */
double VelDiff; /* bandpass selection */
static int first = 0; /* test if first */
/* >>>>>>>> Keplerian radial velocities <<<<<<<<<<<<<<<<<<<< */
axis_1 = Orbit.TrueDistance / (1.0 + 1.0/Binary[Primary].Mq);
if (Flags.elliptic == OFF) {
Orbit.Nu = -0.5 * PI + Phase * (2.0 * PI)/PhaseSteps;
}
FluxOut[Phase].RV[Primary] =
(float) (
(axis_1 / Orbit.TruePeriod)
* 2.0 * PI * sin(Orbit.Inclination)
* (cos(Orbit.Nu + DTOR * Orbit.Omega)
+ Orbit.Excentricity * cos(DTOR * Orbit.Omega))
/ sqrt(1.0 - (Orbit.Excentricity*Orbit.Excentricity))
);
FluxOut[Phase].RV[Secondary] =
(float) ( -(FluxOut[Phase].RV[Primary] / Binary[Primary].Mq) );
/* -------- half amplitudes -----------------------------------*/
Binary[Primary].K =
(axis_1 / Orbit.TruePeriod) * 2.0 * PI * sin(Orbit.Inclination)
/ sqrt(1.0 - (Orbit.Excentricity*Orbit.Excentricity));
Binary[Primary].K = fabs( Binary[Primary].K / 1000. );
Binary[Secondary].K = fabs( Binary[Primary].K / Binary[Primary].Mq );
/* >>>>>>>>>>>>> line profile <<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
if (Flags.lineprofile == ON) {
/* >>>>>>>>>>>> initialize <<<<<<<<<<<<<<<< */
/* ---- the bandpass for the profile -------- */
if (first == 0) {
VelMag = Umag; VelDiff =fabs(1000*WaveL[Primary][Umag]-Orbit.LambdaZero);
for (i = 0; i < NUM_MAG; ++i) {
if (fabs(1000*WaveL[Primary][i] - Orbit.LambdaZero) <= VelDiff) {
VelDiff = fabs(1000*WaveL[Primary][i] - Orbit.LambdaZero);
VelMag = (int) i;
}
}
if (Flags.debug[VERBOSE] == ON)
printf (_(" Line profile computed in %s\n"), Filters[VelMag]);
}
first = 1;
/* ---- the resolution -------- */
DeltaLambda = Orbit.LambdaZero / PROFILE_RES;
/* ---- the conversion factor -------- */
ConstantFactor = Orbit.LambdaZero / (DeltaLambda * SpeedOfLight);
/* ---- the starting wavelength -------- */
Offset = Orbit.LambdaZero - (PROFILE_ARRAY/2) * DeltaLambda;
/* ---- the array for the profile -------- */
for (i = 0; i < PROFILE_ARRAY; ++i) {
Lambda[i] = Offset + i * DeltaLambda;
Profile[Primary][i] = 0.0;
Profile[Secondary][i] = 0.0;
#ifdef HAVE_DISK
Profile[Disk][i] = 0.0; /* Added 2002/04/09, MPR */
#endif
}
/* >>>>>>>>>>>> loop <<<<<<<<<<<<<<<< */
/* Changed for disk component (2->3) 2002/04/09, MPR */
for (Comp = 0; Comp < NUM_COMP; ++Comp) {
SumOfWeights = 0.0;
for (i = 0; i < Binary[Comp].NumElem; ++i) {
if ( Surface[Comp][i].visibility >= FLT_EPSILON ) {
Cosgamma = Surface[Comp][i].CosGamma;
Common = Surface[Comp][i].visibility
* Cosgamma;
OneMinCosgamma = (1.0 - Cosgamma);
/* ---- compute the weights -------- */
if (Flags.limb == 0)
Weight = Surface[Comp][i].f_[VelMag] * Common
*(1.0 - Limb[Comp][VelMag][0]*OneMinCosgamma);
else if (Flags.limb == 1)
Weight = Surface[Comp][i].f_[VelMag] * Common
*(1.0 - Limb[Comp][VelMag][0]*OneMinCosgamma
- Limb[Comp][VelMag][1]*(1.0-SQR(OneMinCosgamma)));
else if (Flags.limb == 2)
Weight = Surface[Comp][i].f_[VelMag] * Common
*(1.0 - Limb[Comp][VelMag][0]*OneMinCosgamma
- Limb[Comp][VelMag][1]*(1.0-sqrt(Cosgamma)));
else if (Flags.limb == 3 && Comp < Disk)
Weight = Surface[Comp][i].f_[VelMag] * Common
*(1.0 - Surface[Comp][i].AreaType * Limb[Comp][VelMag][0]*OneMinCosgamma);
else if (Flags.limb == 3 && Comp >= Disk)
Weight = Surface[Comp][i].f_[VelMag] * Common
*(1.0 - Limb[Comp][VelMag][0]*OneMinCosgamma);
/* ---- Velocity of surface element ----------- */
Velocity = Surface[Comp][i].Velocity + FluxOut[Phase].RV[Comp];
if (Flags.elliptic == OFF && fabs(Velocity) >= FLT_EPSILON)
Velocity = -(Velocity);
/* ---- weighting --------------------------- */
HashLo = (int) (floor( (Velocity*ConstantFactor) ));
HashHi = HashLo + 1;
Lo = fabs((Velocity*ConstantFactor)-HashLo);
Hi = fabs(1.0 - Lo);
if (Lo >= EPS) {
WeightLo = Weight * ( (Hi/Lo) / (1.0 + Hi/Lo) );
WeightHi = Weight - WeightLo;
} else {
WeightLo = Weight;
WeightHi = 0.0;
}
HashLo = (PROFILE_ARRAY/2) + HashLo;
HashHi = (PROFILE_ARRAY/2) + HashHi;
/* ---- put profile in array ------------------- */
if (HashLo >= 0 && HashLo < PROFILE_ARRAY) {
Profile[Comp][HashLo] = Profile[Comp][HashLo] + WeightLo;
SumOfWeights = SumOfWeights + WeightLo;
} else {
WARNING(_("RadVel: Line Profile out of computed Window"));
}
if (HashHi >= 0 && HashHi < PROFILE_ARRAY) {
Profile[Comp][HashHi] = Profile[Comp][HashHi] + WeightHi;
SumOfWeights = SumOfWeights + WeightHi;
} else {
WARNING(_("RadVel: Line Profile out of computed Window"));
}
}
} /* ---- end loop over surface elements ----- */
/* -------- Normalize and Invert the Profile ---- */
if (fabs(SumOfWeights) >= FLT_EPSILON) {
for (i = 0; i < PROFILE_ARRAY; ++i) {
if (fabs(Profile[Comp][i]) >= (FLT_EPSILON) )
Profile[Comp][i] = -fabs(Profile[Comp][i] / SumOfWeights);
}
} else {
for (i = 0; i < PROFILE_ARRAY; ++i) Profile[Comp][i] = 0.0;
}
} /* ------ END loop over components ----- */
/* >>>>>>>>>>>>>>> write output <<<<<<<<<<<<<<<<<<<<<<<<<<<< */
CorrPhase = (FluxOut[Phase].Phase/(PI+PI) - 0.5);
if (CorrPhase <= (-FLT_EPSILON) ) CorrPhase = CorrPhase + 1.0;
sprintf(Out_File, "NightfallProfile_%4.2f.dat", CorrPhase);
file_out = fopen(Out_File, "w");
if (file_out == NULL)
#ifndef _WITH_GTK
nf_error(_(errmsg[3]));
#else
{
if (Flags.interactive == OFF) {
nf_error(_(errmsg[3]));
} else {
make_dialog(_(errmsg[3]));
}
}
#endif
/* ------ file header -------- */
OutfileHeader(file_out);
/* ------- the data ---------- */
fprintf(file_out, "# \n");
fprintf(file_out, _("# PHASE %8.4f\n"), CorrPhase);
fprintf(file_out, "# \n");
fprintf(file_out, _("# Lambda Both Primary Secondary\n") );
for (i = 0; i < PROFILE_ARRAY; ++i) {
temp = Profile[Primary][i] + Profile[Secondary][i];
fprintf(file_out, "%8.4f %11.4g %11.4g %11.4g\n",
Lambda[i], temp,
Profile[Primary][i], Profile[Secondary][i]);
/* UND SOLL DIE DISK AUCH DAZU??? 2002/04/09, MPR */
/* ------- to array ---------- */
ProData[Phase][i] = (float) temp;
}
(void) fclose(file_out);
if (Flags.elliptic == OFF) {
ProData[Phase][PROFILE_ARRAY] =
(float) ( -0.5 + (Orbit.Phase/(2.0*PI)));
} else {
ProData[Phase][PROFILE_ARRAY] =
(float) (-0.5 + ((Orbit.M + PI + Orbit.OmegaZero)/(2.0*PI)));
}
Flags.ProComputed = ON;
}
return;
}
syntax highlighted by Code2HTML, v. 0.9.1