/* 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 #include #include #include #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; }