/* 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 Mark surface areas as spotted @long Computes angular distance to spot centre on great circle, marks surface element as spotted if distance smaller than spot radius. Then temperature is updated according to mean dimfactor (in case of overlapping spots) for a surface element. @param (int Comp) The stellar component @param (int j) The phase index @return (void) @heading Plotting ****************************************************************************/ void MakeSpots(int Comp, int j) { int nspot; /* index of current spot */ long NElem; /* number of surface elements */ int N_Spot = 0; /* # of spots on star */ long i; /* loop counter */ double LatS; /* spot latitude */ double CosLat; /* cosine of spot latitude */ double CosLong; /* cosine of spot longitude */ double SinLat; /* sine of spot latitude */ double EtaS, PhiS; /* eta, phi of spot */ double SinEta; /* sine of eta (spot) */ double CosEta; /* cosine of eta (spot) */ double CosL; /* spheric. trig. stuff */ double DeltaPhi; /* phi(spot) - phi(surf.elem.) */ double SuCosEta, SuSinEta; /* cos/sin eta(surf.elem.) */ double CosRad; /* cosine of spot radius */ double RadSpot; /* spot radius */ double Test = 0.0; /* phase in eccentric orbit */ SurfaceElement *SurfPtr; /* pointer to surface */ double thisDimfactor; /* dimfactor of current spot */ double fratio = 1.0; /* async rotation */ double argA; /* argument for asin */ /* ---------------- initialize -------------------------------------- */ NElem = Binary[Comp].NumElem; if (Comp == Primary) { N_Spot = Flags.Spots1; if (Flags.asynchron1 == ON) fratio = Binary[Primary].Fratio; } else if (Comp == Secondary) { N_Spot = Flags.Spots2; if (Flags.asynchron2 == ON) fratio = Binary[Secondary].Fratio; } SurfPtr = Surface[Comp]; for (i = 0; i < NElem; ++i) { SurfPtr->SpotNum = 0; SurfPtr->SumDimFactor = 0.0; ++SurfPtr; } if (Flags.elliptic == ON) Test = 0.5*PI + (j * (PI+PI)/PhaseSteps) - Orbit.Phase; /* ----------------- loop over spots ---------------------------------- */ for (nspot = 0; nspot < N_Spot; ++nspot) { thisDimfactor = Spot[Comp][nspot].dimfactor; if (Flags.debug[VERBOSE] == ON) printf(_("MakeSpot: Spot %5d on Component %3d \n"), nspot+1, Comp+1); /* ------------ move spot if async rotation ------------------------- */ Spot[Comp][nspot].nowlongitude = Spot[Comp][nspot].longitude - RTOD * Test - (fratio - 1.0) * RTOD * ((PI+PI)/PhaseSteps) * j; if (Comp == Secondary) Spot[Comp][nspot].nowlongitude = Spot[Comp][nspot].nowlongitude + 180; Spot[Comp][nspot].nowlongitude = Spot[Comp][nspot].nowlongitude - floor(Spot[Comp][nspot].nowlongitude/360.) * 360.; /* ---------------- latitude must be +/- 90 deg --------------------- */ /* -----------------longitude is clockwise from positive x-axis ------ */ LatS = DTOR * Spot[Comp][nspot].latitude; /* Must be != NULL, else spot is incorrectly placed */ LatS = ( fabs(LatS) >= DBL_EPSILON) ? LatS : DBL_EPSILON; CosLat = cos(LatS); CosLong = cos(DTOR * Spot[Comp][nspot].nowlongitude); CosEta = CosLat*CosLong; EtaS = acos( CLAMP(CosEta, -1.0, 1.0) ); /* sign is correct */ SinEta = sin(EtaS); /* beware of large roundoff errors */ if (fabs(CosEta-1) <= 0.00001) SinEta = 0.00001; SinLat = sin(LatS); argA = SinLat/SinEta; PhiS = asin( CLAMP(argA, -1.0, 1.0) ); /* sign is not correct */ /* ------------- get proper sign ----------------------------------- */ if (Comp == Primary) { if ( LatS >= DBL_EPSILON && Spot[Comp][nspot].nowlongitude >= (180+DBL_EPSILON) ) PhiS = PI-PhiS; /* PhiS is positive before */ if ( LatS <= (-DBL_EPSILON) && Spot[Comp][nspot].nowlongitude >= (180+DBL_EPSILON) ) PhiS = PI-PhiS; /* PhiS is negative before */ } if (Comp == Secondary) { if ( LatS >= DBL_EPSILON && Spot[Comp][nspot].nowlongitude <= (180-DBL_EPSILON)) PhiS = PI-PhiS; if ( LatS <= (-DBL_EPSILON) && Spot[Comp][nspot].nowlongitude <= ( 180-DBL_EPSILON) ) PhiS = PI-PhiS; } /* ------------- spot size ------------------------------------- */ RadSpot = DTOR * Spot[Comp][nspot].radius; CosRad = cos(RadSpot); /* ------------- Loop over Surface Elements ----------------------- */ SurfPtr = Surface[Comp]; for (i = 0; i < NElem; ++i) { if (fabs(EtaS - SurfPtr->eta) <= (RadSpot-DBL_EPSILON) ) { DeltaPhi = PhiS - SurfPtr->phi; SuCosEta = cos(SurfPtr->eta); SuSinEta = sin(SurfPtr->eta); CosL = CosEta*SuCosEta + SinEta*SuSinEta*cos(DeltaPhi); /* ---------------- If within Spot ------------------------------ */ if (CosRad <= (CosL-DBL_EPSILON)) { ++SurfPtr->SpotNum; SurfPtr->SumDimFactor = SurfPtr->SumDimFactor + thisDimfactor; } } ++SurfPtr; } } /* -------------------------- adjust temperature ------------------------ */ SurfPtr = Surface[Comp]; for (i = 0; i < NElem; ++i) { if (SurfPtr->SpotNum > 0) { SurfPtr->SumDimFactor = SurfPtr->SumDimFactor / SurfPtr->SpotNum; SurfPtr->temp = SurfPtr->SumDimFactor * SurfPtr->temp; } else { SurfPtr->SumDimFactor = 1.0; } ++SurfPtr; } return; }