/* 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 #include "Light.h" /* 2002-05-03 rwichmann HIGH_PRECISION code */ #ifdef HIGH_PRECISION /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.2 @short Numerically differentiate dR/deta @param (int) Comp The stellar component @return (double) dR/deta @heading Star Setup ****************************************************************************/ double NumDiffR (double deltaIn, double etaIn, double sinphi, double RadiusBound, double RXCrit, int * testerr) { double val1, val0; double eps = pow(DBL_EPSILON, (1.0/3.0)); double delta; double nu_store; double lambda_store; delta = eps * deltaIn; delta = (delta < eps) ? eps : delta; nu_store = nu; lambda_store = lambda; nu = sinphi * sin(etaIn-delta); lambda = cos(etaIn-delta); val0 = RootFind(RocheSurface, RadiusBound, RXCrit, 100.0 * DBL_EPSILON, "NumDiffR", testerr); if (*testerr == 1) return(0.0); nu = sinphi * sin(etaIn+delta); lambda = cos(etaIn+delta); val1 = RootFind(RocheSurface, RadiusBound, RXCrit, 100.0 * DBL_EPSILON, "NumDiffR", testerr); if (*testerr == 1) return(0.0); nu = nu_store; lambda = lambda_store; return ((val1-val0)/(2.0*delta)); } /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.2 @short Numerically differentiate dR/dphi @param (int) Comp The stellar component @return (double) dR/dphi @heading Star Setup ****************************************************************************/ double NumDiffRPhi (double deltaIn, double phiIn, double sineta, double RadiusBound, double RXCrit, int * testerr) { double val1, val0; double eps = pow(DBL_EPSILON, (1.0/3.0)); double delta; double nu_store; delta = eps * deltaIn; delta = (delta < eps) ? eps : delta; nu_store = nu; nu = sineta * sin(phiIn-delta); val0 = RootFind(RocheSurface, RadiusBound, RXCrit, 100.0 * DBL_EPSILON, "NumDiffRPhi", testerr); if (*testerr == 1) return(0.0); nu = sineta * sin(phiIn+delta); val1 = RootFind(RocheSurface, RadiusBound, RXCrit, 100.0 * DBL_EPSILON, "NumDiffRPhi", testerr); if (*testerr == 1) return(0.0); nu = nu_store; return ((val1-val0)/(2.0*delta)); } /* #ifdef HIGH_PRECISION */ #endif /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.2 @short Divide stellar surface into grid @long First, set up the grid, then compute the temperature distribution on the grid, including spots and reflection @param (int) Comp The stellar component @return (int) The exit status @heading Star Setup ****************************************************************************/ int DivideSurface(int Comp) { register long i,j; /* loop variables */ long n_phi; /* # of surface steps in phi */ long n_eta; /* # of surface steps in eta */ long NumElem; /* current surface element */ long NextPtr; /* next surface element */ int testerr = 0; /* exit status of RootFind */ double eta; /* eta of surface element */ double phi; /* phi of surface element */ double sineta, sinphi; /* sin,cos of phi */ double coseta, cosphi; /* sin,cos of eta */ double rad1; /* radius of surface element */ double sqr_rad1; /* square of radius */ double delta_eta; /* step in eta */ double delta_phi; /* step in phi */ double rho; /* rho of surface element */ /* 2002-05-03 rwichmann HIGH_PRECISION code */ #ifdef HIGH_PRECISION double rdiff; /* dR/deta */ double rdiff2; /* dR/dphi */ #endif double cub_rho; /* cube of rho */ double C_eta, C_phi, C_rad1; /* surface gravity vector */ /* in spherical coordinates */ double C_x, C_y, C_z; /* surface gravity vector */ /* in cartesian coordinates */ double grav; /* surface gravity */ double lvec, mvec, nvec; /* surface normal vector */ double RXCrit; /* Critical radius */ double DSurf; /* surface of element */ double MeanGrav = 0.0; /* mean gravity */ double Surf = 0.0; /* total surface */ double MeanDarkGrav = 0.0; /* contribution to grav. darkening */ double RadiusBound = 0.0; /* inner bound for RootFind */ double T1; /* temperature of this star */ double T2; /* temperature of the other star */ double Tr; /* temperature ratio^4 */ double etacorr_over, eta_corr, deltaeta_corr; /* Throat adjustment */ double eta_next; /* Pointer setting */ long n_next, n_last, n_comp; /* Pointer setting */ /* for mean temperature calculation */ double Area; /* area of element */ double Weight=0.0; /* total area */ double TempMean=0.0; /* mean temperature */ SurfaceElement *SurfPtr; /* pointer to surface */ BinaryComponent *BinPtr; /* pointer to binary */ /* ---------- initialize ---------------------------------------------- */ SurfPtr = Surface[Comp]; BinPtr = &Binary[Comp]; Mq = Binary[Comp].Mq; if (Comp == Primary && Flags.asynchron1 == ON) F = Binary[Comp].Fratio; else if (Comp == Secondary && Flags.asynchron2 == ON) F = Binary[Comp].Fratio; else F = 1.0; RochePot = Binary[Comp].RochePot; RadiusBound = 0.8*Binary[Comp].Radius; etacorr_over = Binary[Comp].LimEta; /* ---------- Maximum Radius to search for Root ---------------- */ if (Flags.fill == OFF) { RXCrit = BinPtr->RXCrit; } else { RXCrit = BinPtr->LimRad + (0.07*SQR(Binary[Primary].Mq)) + 0.02*Binary[Primary].Mq; if ( (Binary[Primary].Mq - 1.0) >= FLT_EPSILON) RXCrit = BinPtr->LimRad + (0.07*SQR(Binary[Secondary].Mq)) + 0.02*Binary[Secondary].Mq; if ( (Binary[Primary].RocheFill - 1.05) <= FLT_EPSILON && (Binary[Primary].Mq - 0.2) >= FLT_EPSILON && (Binary[Primary].Mq - 5.0) <= FLT_EPSILON ) RXCrit = BinPtr->RXCrit + Binary[Primary].RocheFill - 1.; else if ( (Binary[Primary].RocheFill - 1.02) <= FLT_EPSILON ) RXCrit = BinPtr->RXCrit + Binary[Primary].RocheFill - 1.; /* critical angle w/r to z */ etacorr_over = BinPtr->RZatL1 / sqrt( BinPtr->RZatL1 * BinPtr->RZatL1 + BinPtr->RLag1 * BinPtr->RLag1); etacorr_over = asin( CLAMP(etacorr_over, -1.0, 1.0) ); } BinPtr->Surface = 0.0; BinPtr->Gravity = 0.0; if (Comp == 0) { T1 = BinPtr->Temperature; T2 = Binary[Comp+1].Temperature; } else { T1 = BinPtr->Temperature; T2 = Binary[Comp-1].Temperature; } Tr = T2/T1; Tr = Tr * Tr * Tr * Tr; n_eta = StepsPerPi; n_phi = 0; if (Flags.fill == OFF) BinPtr->LimEta = 0.0; delta_eta = (PI - BinPtr->LimEta)/n_eta; eta = PI-(0.5*delta_eta); /* start value, point opposite to companion */ if (Flags.fill == ON) { /* Throat adjustment */ etacorr_over = (BinPtr->LimEta - etacorr_over) + 3.0 * delta_eta; } NumElem = 0; /* ------------------- loop over eta ------------------------------ */ BinPtr->N_Rings = n_eta; for (i = 0; i < n_eta; ++i) { eta_corr = eta; /* Throat adjustment */ deltaeta_corr = delta_eta; phi = 0.0; sineta = sin(eta); coseta = cos(eta); lambda = coseta; /* fix n_phi for GL animation, because the vertex list needs a */ /* fixed number vertices per phi loop */ #ifdef _WITH_OPENGL if (Flags.animate == ON && Flags.use_opengl == ON) { /* n_phi fixed to a power of 2 for GL animation */ n_phi = 32; } else { /* n_phi proportional to sin(eta) */ n_phi = (long) (10 + 1 + (int)(fabs(2 * n_eta * sineta))); } #else /* n_phi proportional to sin(eta) */ n_phi = (long) (10 + 1 + (int)(fabs(2 * n_eta * sineta))); #endif delta_phi = (PI+PI)/n_phi; /* check if enough memory */ if ((NumElem+n_phi) >= (int) MaximumElements) { #ifdef _WITH_GTK if (Flags.interactive == OFF) nf_error(_(errmsg[10])); else { if (Flags.WantFit == OFF && Flags.WantMap == OFF) make_dialog(_(errmsg[10])); return(8); } #else nf_error(_(errmsg[10])); #endif } BinPtr->N_PhiStep[i] = n_phi; /* store number of phi intervals */ /* --------- Loop over phi ---------------------------------------- */ for (j = 0; j < n_phi; ++j) { sinphi = sin(phi); cosphi = cos(phi); /* Throat adjustment for last three rings */ /* The throat is wider in x-y plane, thus we need to adjust eta, */ /* otherwise there would be an increasing gap towards higher */ /* latitudes. */ if ((Flags.fill == ON) && (i > (n_eta-4))) { if (i == n_eta-3) eta_corr = eta - ((etacorr_over/6.0) -(delta_eta/2.0))*fabs(sinphi); else if (i == n_eta-2) eta_corr = eta - ((etacorr_over/2.0) -(1.5*delta_eta))*fabs(sinphi); else if (i == n_eta-1) eta_corr = eta - ((5.0*etacorr_over/6.0) -(5.0*delta_eta/2.0))*fabs(sinphi); sineta = sin(eta_corr); coseta = cos(eta_corr); lambda = coseta; deltaeta_corr = delta_eta+((etacorr_over/3.0)-delta_eta)*fabs(sinphi); } mu = sineta*cosphi; nu = sineta*sinphi; rad1 = RootFind(RocheSurface, RadiusBound, RXCrit, 100.0 * DBL_EPSILON, /* was 1.0e-8, */ "DivideSurface1", &testerr); if (testerr == 1) return(8); sqr_rad1 = rad1*rad1; rho = 1.0/sqrt(1.0+(sqr_rad1)-2.0*rad1*lambda); cub_rho = rho*rho*rho; /* old code C_rad1 = -1.0/(sqr_rad1) + Mq*(cub_rho*(lambda-rad1)-lambda) +F*F*(Mq+1.0)*(sqr_rad1)*(1.0- nu*nu); */ /* 2002-05-03 rwichmann error in expression fixed */ C_rad1 = -1.0/(sqr_rad1) + Mq*(cub_rho*(lambda-rad1)-lambda) +F*F*(Mq+1.0)*rad1*(1.0- nu*nu); C_eta = sineta *(Mq*(1.0-cub_rho)-(Mq+1.0)*rad1*F*F *lambda*sinphi*sinphi); C_phi = -(Mq + 1.0)*F*F*rad1*nu*cosphi; C_x = C_rad1*lambda - C_eta*sineta; C_y = C_rad1*mu + C_eta*lambda*cosphi - C_phi*sinphi; C_z = C_rad1*nu + C_eta*lambda*sinphi + C_phi*cosphi; grav = sqrt( C_x*C_x + C_y*C_y + C_z*C_z ); lvec = -C_x/grav; mvec = -C_y/grav; nvec = -C_z/grav; DSurf = (sqr_rad1) * sineta * deltaeta_corr * delta_phi; /* 2002-05-03 rwichmann HIGH_PRECISION code */ #ifdef HIGH_PRECISION /* * dA = sin(eta) deta dphi * ( R*R e(r) - R dR/deta e(eta) * - R/sin(eta) dR/dphi e(phi) ) */ /* the dR/deta term */ rdiff = NumDiffR (deltaeta_corr, eta_corr, sinphi, RadiusBound, RXCrit, &testerr); if (testerr == 1) return (8); rdiff = rad1 * rdiff * sineta * deltaeta_corr * delta_phi; /* the dR/dphi term */ rdiff2 = NumDiffRPhi (delta_phi, phi, sineta, RadiusBound, RXCrit, &testerr); if (testerr == 1) return (8); rdiff2 = rad1 * rdiff2 * deltaeta_corr * delta_phi; /* absolute value */ DSurf = sqrt(DSurf*DSurf + rdiff*rdiff + rdiff2*rdiff2); #else /* * cos (e) = lambda*lvec+mu*mvec+nu*nvec * = angle to surface normal */ DSurf = DSurf/(lambda*lvec+mu*mvec+nu*nvec); #endif SurfPtr->eta = (float) eta_corr; SurfPtr->phi = (float) phi; SurfPtr->rad = (float) rad1; SurfPtr->lambda = (float) (lambda * rad1); SurfPtr->mu = (float) (mu * rad1); SurfPtr->nu = (float) (nu * rad1); SurfPtr->rho = (float) rho; SurfPtr->grav = (float) grav; SurfPtr->l = (float) lvec; SurfPtr->m = (float) mvec; SurfPtr->n = (float) nvec; SurfPtr->area = (float) DSurf; SurfPtr->MyNum = NumElem; SurfPtr->ring = i; Surf = Surf + DSurf; MeanGrav = MeanGrav + (grav*DSurf); MeanDarkGrav = MeanDarkGrav + exp(BinPtr->GravDarkCoeff * log(grav))*DSurf; /* --------------- INITIALIZE SURFACE TEMP ----------------- */ /* If we want to consider limb brightening for illuminated side * of star, we need to mark illuminated/non-illuminated surface * elements. Let the default be non-illuminated, and mark illuminated * surface elements in LightReflect(). */ SurfPtr->AreaType = BACKSIDE; SurfPtr->temp = T1; phi = phi + delta_phi; /* increment phi */ ++NumElem; /* increment surface element count */ ++SurfPtr; } /* ----------- end loop over phi --------------------------- */ eta = eta - delta_eta; } /* -------------- end loop over eta ------------------------------- */ BinPtr->NumElem = NumElem; BinPtr->NumPlus = NumElem; if (Flags.debug[VERBOSE] == ON) { printf("\n"); printf(_(" Divide Star %3d : %6ld Surface elements\n"), Comp+1, BinPtr->NumElem); } BinPtr->Surface = Surf; BinPtr->Gravity = MeanGrav/Surf; BinPtr->DarkeningGrav = MeanDarkGrav/Surf; /* ------------------- SET next/last_ptr ---------------------------- */ /* we need this only once */ if (Flags.first_pass == ON) { SurfPtr = Surface[Comp]; n_last = 0; n_next = 0; NumElem = 0; for (i = 0; i < n_eta; ++i) { eta_next = eta - delta_eta; /* pointer setting */ if (i > 0) n_last = n_phi; /* pointer setting */ n_phi = BinPtr->N_PhiStep[i]; if ( i < n_eta-1) n_next = BinPtr->N_PhiStep[i+1]; for (j = 0; j < n_phi; ++j) { /* ----------- the first ring ---------------------------- */ if (i == 0) { SurfPtr->last_ptr = NULL; if (n_next == n_phi) { SurfPtr->next_ptr = &Surface[Comp][NumElem+n_next]; } else { n_comp = ROUND(((float)(n_next)/(float)(n_phi))*(float)(j)); n_comp = n_phi - j + n_comp; SurfPtr->next_ptr = &Surface[Comp][NumElem+n_comp]; } } /* ----------- somewhere in the middle -------------------- */ else if ( (i > 0) && (i < (n_eta-1))) { if (n_next == n_phi) { NextPtr = NumElem+n_next; if (NextPtr >= BinPtr->NumElem) NextPtr = BinPtr->NumElem-1; SurfPtr->next_ptr = &Surface[Comp][NextPtr]; } else { n_comp = ROUND(((float)(n_next)/(float)(n_phi))*(float)(j)); n_comp = n_phi - j + n_comp; NextPtr = NumElem+n_comp; if (NextPtr >= BinPtr->NumElem) NextPtr = BinPtr->NumElem-1; SurfPtr->next_ptr = &Surface[Comp][NextPtr]; } if (n_last == n_phi) { SurfPtr->last_ptr = &Surface[Comp][NumElem-n_next]; } else { n_comp = ROUND(((float)(n_last)/(float)(n_phi))*(float)(j)); n_comp = (n_last - n_comp) + j; SurfPtr->last_ptr = &Surface[Comp][NumElem-n_comp]; } } /* ----------- the last ring ---------------------------- */ /*@ignore@*/ /* Storage SurfPtr->last_ptr is kept in one path, but live in */ /* another. */ else if (i == (n_eta-1)) { if (n_last == n_phi) { NextPtr = NumElem-n_phi; if (NextPtr >= BinPtr->NumElem) NextPtr = BinPtr->NumElem-1; SurfPtr->last_ptr = &Surface[Comp][NextPtr]; } else { n_comp = ROUND(((float)(n_last)/(float)(n_phi))*(float)(j)); n_comp = (n_last - n_comp) + j; NextPtr = NumElem-n_comp -1; if (NextPtr >= BinPtr->NumElem) NextPtr = BinPtr->NumElem-1; SurfPtr->last_ptr = &Surface[Comp][NextPtr]; } if (Flags.fill == OFF) { n_comp = (j + (n_phi/2)) % n_phi; n_comp = - j + n_comp; NextPtr = NumElem+n_comp; SurfPtr->next_ptr = &Surface[Comp][NextPtr]; } else { if (Comp == Primary) { SurfPtr->next_ptr = &Surface[Secondary][NextPtr]; } else { NextPtr = NumElem; SurfPtr->next_ptr = &Surface[Primary][NextPtr]; } } } if (NumElem != 0) { if (j != 0) { SurfPtr->phi_ptr = &Surface[Comp][NumElem-1]; } else { SurfPtr->phi_ptr = &Surface[Comp][NumElem+n_phi-1]; } if (j != (n_phi-1)) { SurfPtr->phi1_ptr = &Surface[Comp][NumElem+1]; } else { SurfPtr->phi1_ptr = &Surface[Comp][NumElem-(n_phi-1)]; } } else { SurfPtr->phi1_ptr = &Surface[Comp][1]; SurfPtr->phi_ptr = &Surface[Comp][1]; } ++NumElem; ++SurfPtr; } /*@end@*/ } } /* -------------- END SET next/last_ptr ----------------------------- */ /**********************************************************************/ /* */ /* Temperature Distribution */ /* */ /**********************************************************************/ /* --------------- Gravitational Darkening ------------------------- */ /* T propto g**GravDarkCoeff */ SurfPtr = Surface[Comp]; for (i = 0; i < BinPtr->NumElem; ++i) { SurfPtr->temp = BinPtr->Temperature * (exp(BinPtr->GravDarkCoeff * log(SurfPtr->grav)) / BinPtr->DarkeningGrav ); SurfPtr->orig_temp = SurfPtr->temp; ++SurfPtr; } /* ------------------------ Spots -------------------- */ /* make spots or initialize the dimfactor to 1.0 */ if (Flags.first_pass == ON) { if (Comp == Primary) { if (Flags.Spots1 > OFF) { MakeSpots(Primary, 0); } else { SurfPtr = Surface[Primary]; for (j = 0; j < Binary[Primary].NumElem; ++j) { SurfPtr->SumDimFactor = 1.0; ++SurfPtr; } } } if (Comp == Secondary) { if (Flags.Spots2 > OFF) { MakeSpots(Secondary, 0); } else { SurfPtr = Surface[Secondary]; for (j = 0; j < Binary[Secondary].NumElem; ++j) { SurfPtr->SumDimFactor = 1.0; ++SurfPtr; } } } } else { /* move spots if asynchroneous rotation or elliptic orbit */ if (Comp == Primary) { if ( (Flags.Spots1 > OFF) && ( Flags.asynchron1 == ON || Flags.elliptic == ON) ) MakeSpots(Primary, Orbit.PhaseIndex); } if (Comp == Secondary) { if ( (Flags.Spots2 > OFF) && ( Flags.asynchron2 == ON || Flags.elliptic == ON) ) MakeSpots(Secondary, Orbit.PhaseIndex); } } /* --------------------------------- Reflection -------------------- */ if ((Flags.reflect > OFF) && (Comp == Secondary)) { LightReflect(); } if (Flags.reflect == OFF) SimpleReflect(Comp); /* ---------------------------- mean temperature ------------------- */ Weight = 0.0; TempMean = 0.0; SurfPtr = Surface[Comp]; for (i = 0; i < BinPtr->NumElem; ++i) { Area = SurfPtr->area; TempMean = TempMean + SurfPtr->temp * Area; Weight = Weight + Area; /*------- Sollte hier nicht irgendwo ein SurfPtr++ stehen????? */ /* eigentlich schon ... RW Tue Sep 11 12:54:20 CEST 2001 */ ++SurfPtr; } BinPtr->TempMean = TempMean/Weight; /**********************************************************************/ /* */ /* Flux intensity computation */ /* */ /**********************************************************************/ if (Comp == Secondary) { if (Flags.blackbody == ON) { BlackbodyFlux(Primary); BlackbodyFlux(Secondary); } else { ModelFlux(Primary); ModelFlux(Secondary); } } /* ---------------- output --------------------- */ if (Flags.debug[VERBOSE] == ON) { /* printf("\n Star: %5d\n", (Comp+1)); */ printf(_(" Surface Area, Mean Gravity (dimensionless): %8.6f %8.6f\n"), BinPtr->Surface, BinPtr->Gravity); printf(_(" Temperature is %8.6f\n"), BinPtr->Temperature); printf("\n"); printf(_(" Temperature Distribution corrected for: \n")); printf(_(" Reflection and Gravity Darkening\n")); printf(_(" Coefficient for Gravity Darkening was: %6.4f\n"), BinPtr->GravDarkCoeff); } /* ---------------- statistics --------------------- */ if (Comp == Secondary && Flags.debug[SURFACE] == ON) { for (i = 0; i < 2; ++i) { int ring; int NumElem = Binary[i].NumElem; double r_temp[STEPS_PER_PI] = { 0.0 }; double r_flux[STEPS_PER_PI] = { 0.0 }; double r_refl[STEPS_PER_PI] = { 0.0 }; int r_n[STEPS_PER_PI] = { 0 }; for (j = 0; j < NumElem; ++j) { ring = Surface[i][j].ring; r_temp[ring] += Surface[i][j].orig_temp; r_refl[ring] += Surface[i][j].temp; r_flux[ring] += Surface[i][j].f_[Rmag]; ++r_n[ring]; } printf("\n"); printf("%02ld Temp Trefl Ratio Flux\n", i+1); printf("------------------------------------------------------\n"); for (j = 0; j < Binary[i].N_Rings; ++j) { printf("%02ld %10.4g %10.4g %10.4g %10.4g %03d\n", j, r_temp[j]/r_n[j], r_refl[j]/r_n[j], r_refl[j]/r_temp[j], r_flux[j]/r_n[j], r_n[j]); } printf("\n"); } } return(0); } #if 0 /* Can we rename this to a more intuitive name light rotatex or what ever; I don't know what "tile" in this context means; Why using floats and not double values ? MK */ /* This is a no-op if the tilt angle is zero */ static void TileDisk1(double TiltAngle,SurfaceElement *SurfPtr) { float lambda_strich; float nu_strich; lambda_strich = (float) SurfPtr -> lambda*cos(TiltAngle)- SurfPtr -> nu*sin(TiltAngle); nu_strich = (float) SurfPtr -> lambda*sin(TiltAngle)+ SurfPtr -> nu*cos(TiltAngle); SurfPtr->lambda = lambda_strich; SurfPtr->nu = nu_strich; } /* This is a no-op if the tilt angle is zero */ static void TileDisk2(double TiltAngle,SurfaceElement *SurfPtr) { float l_strich; float n_strich; l_strich = (float) SurfPtr -> l*cos(TiltAngle)+ SurfPtr -> n*sin(TiltAngle); n_strich = (float) SurfPtr -> l*sin(TiltAngle)- SurfPtr -> n*cos(TiltAngle); SurfPtr -> l = l_strich; SurfPtr -> n = n_strich; } /* This is a no-op if the tilt angle is zero */ static void TileDisk3(double TiltAngle,SurfaceElement *SurfPtr) { float lambda_strich; float mu_strich; lambda_strich = (float) SurfPtr -> lambda*cos(TiltAngle)- SurfPtr -> mu*sin(TiltAngle); mu_strich = (float) SurfPtr -> lambda*sin(TiltAngle)- SurfPtr -> mu*cos(TiltAngle); SurfPtr -> lambda = lambda_strich; SurfPtr -> mu = mu_strich; } static void CalculateLayerVecE1(SurfaceElement *SurfPtr) { double vec_en[3]; /* Normalenvektor der Ebene */ vec_en[0] = SurfPtr->l; vec_en[1] = SurfPtr->m; vec_en[2] = SurfPtr->n; if (vec_en[0] != 0) { SurfPtr->VecE1[0] = -((vec_en[1]+vec_en[2])/vec_en[0]); SurfPtr->VecE1[1] = 1.0; SurfPtr->VecE1[2] = 1.0; } else if (vec_en[1] != 0) { SurfPtr->VecE1[0] = 1.0; SurfPtr->VecE1[1] = -((vec_en[0]+vec_en[2])/vec_en[1]); SurfPtr->VecE1[2] = 1.0; } else { SurfPtr->VecE1[0] = 1.0; SurfPtr->VecE1[1] = 1.0; SurfPtr->VecE1[2] = -((vec_en[0]+vec_en[1])/vec_en[2]); } } static void CalculateLayerVecE2(SurfaceElement *SurfPtr) { double vec_en[3]; /* Normalenvektor der Ebene */ double vec_e1[3]; /* Vektor Ebene 1 */ vec_en[0] = SurfPtr->l; vec_en[1] = SurfPtr->m; vec_en[2] = SurfPtr->n; vec_e1[0] = SurfPtr->VecE1[0]; vec_e1[1] = SurfPtr->VecE1[1]; vec_e1[2] = SurfPtr->VecE1[2]; SurfPtr->VecE2[0] = vec_en[1]*vec_e1[2]-vec_en[2]*vec_e1[1]; SurfPtr->VecE2[1] = vec_en[2]*vec_e1[0]-vec_en[0]*vec_e1[2]; SurfPtr->VecE2[2] = vec_en[0]*vec_e1[1]-vec_en[1]*vec_e1[0]; } #endif #ifdef HAVE_DISK /* This macro yields correct results only if a > 0 and b in [0, 2*PI). */ #define DIST_CIRCLE(a, b) fabs(PI - fabs(PI - fmod((a) - (b), (2.0 * PI) ) )) /**************************************************************************** @author Patrick Risse @short Divide disk surface into grid @long First, set up the grid, then compute the temperature distribution on the grid, including spots and reflection @param (int) Comp The stellar component @return (int) The exit status @heading Disk Setup ****************************************************************************/ int DivideDisk(int Comp) { register long i,j; /* loop variables */ int n_phi = 32.0; /* surface steps in phi */ long NumElem; /* current surface element */ /*long NextPtr;*/ /* next surface element */ /*int testerr = 0;*/ /* exit status of RootFind */ double phi; /* phi of surface element */ double RadiusBound = 0.0; /* inner bound for RootFind */ double T1; /* temperature of this star */ double R1; /* Roche fill of this star */ double Rout,Rin; /* Outer/inner Radius of the Disk */ double etacorr_over; /* Throat adjustment */ /*long n_next, n_last, n_comp;*/ /* Pointer setting */ /* for mean temperature calculation */ double Area; /* area of element */ double Weight=0.0; /* total area */ double TempMean=0.0; /* mean temperature */ double Phase; /* Phase of the System. Is needed for recalculating the Disk Position*/ SurfaceElement *SurfPtr; /* pointer to surface */ BinaryComponent *BinPtr; /* pointer to binary */ /* ---------------- Variables added for the disk by P.R. ---------------*/ double Thick_out; /* Thickness of the disk on the outside */ double Thick_in; /* Thickness of the disk on the inside */ double diskExp = 1.0; /* Exponent for r dependency of scale height */ double diskProp = 1.0; /* Proportionality factor for disk height */ double diskAdd = 1.0; /* Additive constant for disk height */ double d_s1d; /* Area of one element on the top/ bottom side of the disk */ double d_s2d_out; /* Area */ double d_s2d_in; double P_0; double DSurf; /* surface of element */ double DSurfCor; /* surface of element, sin included */ double Surf = 0.0; /* total surface */ double ring_height; double r_i_old; double r_i_new; double h_i_old; double h_i_new; double rise_angle; double r_inner_sqr; double r_outer_sqr; double dphi = 0.0; double r_j=0; double alpha; /*Angles in Degree*/ double TiltAngle; double offset_in=0; double offset_out=0; int n_r; /* Number of Rings on the top/bottom of the disk */ /*long TopElem,OutElem,BottomElem,InElem;*/ int NumOfTopElements,NumOfOuterElements ; int NumOfBottomElements,NumOfInnerElements; int NumOfTotalElements; double warpfactor[42]; /* Inclination of each ring for the warp*/ float rot[42]; int warp_no; /* Entry_No of warpfactor */ int tilt_flag = 0; double t_c; double t_temp; int t_flag = Flags.tdisk; double sin_phi, cos_phi, sin_alpha, cos_alpha; double hs_long = fmod(DTOR * Binary[Disk].longitudeHS, 2.0 * PI); /* longitude of hot spot */ double hs_ext = ((Binary[Disk].extentHS / 2.0) * DTOR); /* extent in phi */ double hs_depth = Binary[Disk].depthHS; /* depth in R */ double hs_temp = Binary[Disk].tempHS; /* temperature */ double hs_phi; /* temporary variable */ /* If hs_long is negative, map it into [0, 2.0 * PI), otherwise * the DIST_CIRCLE algorithm does not work. */ hs_long = ((hs_long < 0) ? (2.0 * PI + hs_long) : hs_long); n_r = floor(StepsPerPi/1.5); if (t_flag == 0) /* simple */ diskExp = 1.0; else if (t_flag == 1) /* isothermal */ diskExp = 1.5; else if (t_flag == 2) /* reprocessing */ diskExp = (9.0/8.0); /*BinPtr->N_PhiStep[i] = n_phi;*/ /* store number of phi intervals */ /* ---------- initialize ---------------------------------------------- */ SurfPtr = Surface[Comp]; BinPtr = &Binary[Comp]; Mq = Binary[Comp].Mq; if (Comp == Primary && Flags.asynchron1 == ON) F = Binary[Comp].Fratio; else if (Comp == Secondary && Flags.asynchron2 == ON) F = Binary[Comp].Fratio; else F = 1.0; /* RLag1 = Binary[Comp].RLag1; * This Value is not used. Why is it defined? P.R. 19.6.2001 */ RochePot = Binary[Comp].RochePot; RadiusBound = 0.8*Binary[Comp].Radius; etacorr_over = Binary[Comp].LimEta; Phase = 0.0; /* ((2*PI/PhaseSteps) * index) */; BinPtr->RocheFill = Binary[Secondary].RocheFill; BinPtr->RCrit = Binary[Secondary].RCrit; BinPtr->RLag1 = Binary[Secondary].RLag1; /* RLag1 is not calculated for the Disk so i take it from the * Secondary Star. There must be a better possibility to define * the RLag1 for the Disk. ToDo: Where is RLag1 defined? * P.R. 19.06.2001 */ Rout = BinPtr->RLag1*BinPtr->Rout; Rin = BinPtr->RLag1*BinPtr->Rin; /*BinPtr->RocheFill; */ BinPtr->Surface = 0.0; BinPtr->Gravity = 0.0; T1 = BinPtr->Temperature; R1 = BinPtr->Radius; NumElem = 0; BinPtr->Segments = n_phi; /* We need to define the proportionality factor such that the */ /* disk temperature matches at Rin: T1 = t_c * pow(Rin, -0.75) */ t_c = (T1 / pow(Rin, -0.75)); /* We need to define the proportionality factor such that the disk */ /* height matches at Rin: Thick = diskProp * pow(Rin, diskExp) */ if (t_flag == 0) { diskAdd = BinPtr->Thick; diskProp = BinPtr->HR; } else { diskAdd = 0.0; diskProp = (BinPtr->Thick / pow(Rin, diskExp)); } /* damn is this a hard job ! It's difficult to set the number of */ /* phi steps to n_phi, this finally took 1 year to be done MK */ /* BinPtr->N_PhiStep[0] = n_phi; */ /* This Function is not tested. There are some bugs inside */ /* I will fix it later MPR 12.09.2001 */ if (Flags.warp == ON ) { /* Create warpfactor */ printf("Creat standard warp\n"); for (i=0; i< n_r-1; i++) { if (i <= 15) warpfactor[i]=(15.0/360)*2*PI; if (i > 15) warpfactor[i]=((15.0 -((i-15)*0.75))/360)*2*PI; if (warpfactor[i] <= 0) warpfactor[i]=0.0; if (i >= 1) rot[i]=((2.5/360)*2*PI)*i; } } /* ------------------- Start building Disk grid --------------------- */ /* --------------------- General Calculations ------------------------ */ /* HR := Height/Radius * phi := angle in xy-plane * n_phi := no. of sectors * n_r := no. of rings */ /* Thickness of the Disk on the outer side */ Thick_out = diskAdd + diskProp * pow(Rout, diskExp); /* Thickness of the Disk on the inner side */ Thick_in = diskAdd + diskProp * pow(Rin, diskExp); /* Average area of one element */ d_s1d = (PI*((Rout*Rout)-(Rin*Rin)))/(n_r*n_phi); /* Area of one element on the outer side */ d_s2d_out = (2.0*PI*Rout*Thick_out)/n_phi; /* Area of one element on the inner side */ d_s2d_in = (2.0*PI*Rin*Thick_in)/n_phi; /* Number of rings rounded (to make elements approx. equal sized) */ BinPtr->RingsOut = floor((d_s2d_out/d_s1d)+0.5); /* Number of rings rounded (to make elements approx. equal sized) */ BinPtr->RingsIn = floor((d_s2d_in/d_s1d)+0.5); P_0 = (Rout*Rout-Rin*Rin)/n_r; r_i_old = Rin; TiltAngle = (BinPtr->Tilt/360)*2*PI; if (fabs(TiltAngle) > FLT_EPSILON) tilt_flag = 1; if (Flags.debug[GEO]) { printf("\nDisk Geometry\n----------------------------------------\n"); printf("Sectors : %d \n", n_phi); printf("Rings : %d \n", n_r); printf("RingsOut: %d \n", BinPtr->RingsOut); printf("ThickOut: %f \n", Thick_out); printf("RingsIn : %d \n", BinPtr->RingsIn); printf("ThickIn : %f \n", Thick_in); } /* ----------------- Calculate number of elements for the Disk------ */ NumOfTopElements = n_phi * n_r; NumOfOuterElements = BinPtr->RingsOut * n_phi; NumOfBottomElements = n_phi * n_r; NumOfInnerElements = BinPtr->RingsIn * n_phi; NumOfTotalElements = NumOfTopElements + NumOfOuterElements + NumOfBottomElements + NumOfInnerElements; if (Flags.debug[GEO]) { printf("Elements: %d\n", NumOfTotalElements); } /**********************************************************************/ /* */ /* ---------------------- Calculate top of disk -------------------- */ /* */ /**********************************************************************/ r_inner_sqr = (r_i_old*r_i_old); /* Area of an element. This can be lifted outside the loop, since * the progression of radii ensures equal areas for each element. */ DSurf = (PI*P_0)/n_phi; /* ----------------- loop over rings (n_r) ---------------------- */ for (i = 0; i < n_r; i++) { phi = 0.0; /* check if enough memory */ if ((NumElem+n_phi) >= (int) MaximumElements) { if (Flags.debug[GEO]) printf("Top of disk: %ld >= %d\n", (NumElem+n_phi), (int) MaximumElements); #ifdef _WITH_GTK if (Flags.interactive == OFF) nf_error(_(errmsg[10])); else { if (Flags.WantFit == OFF && Flags.WantMap == OFF) make_dialog(_(errmsg[10])); return(8); } #else nf_error(_(errmsg[10])); #endif } /* Radius of next circle. This progression ensures that all rings * have equal surface if the rings are flat. */ r_i_new = sqrt(r_inner_sqr + ((i+1)*P_0)); /* Center of element */ r_j = (r_i_new+r_i_old)/2.0; /* Height of this ring */ ring_height = ( (diskAdd + diskProp * pow(r_j, diskExp)) / 2.0); /* angle of rise for this ring */ h_i_old = ( (diskAdd + diskProp * pow(r_i_old, diskExp)) / 2.0); h_i_new = ( (diskAdd + diskProp * pow(r_i_new, diskExp)) / 2.0); rise_angle = atan ( fabs(h_i_new - h_i_old) / fabs(r_i_new - r_i_old) ); /* Corrected surface area */ DSurfCor = DSurf / cos ( rise_angle ); /* shift rings in phi by random amount to reduce numerical artifacts */ #ifdef _WITH_OPENGL if (Flags.use_opengl == OFF) #endif dphi = (2.0*PI * (rand() / (RAND_MAX + 1.0))); /* Angle of rise from the middle to the edge of the Accretion Disk */ if (t_flag == 0) alpha = 0.5 * diskProp; else alpha = rise_angle /* 0.5 * diskProp * diskExp * pow(r_j, (diskExp - 1.0)) */; sin_alpha = sin(alpha); cos_alpha = cos(alpha); if (t_flag < 2) /* isothermal */ t_temp = T1; else t_temp = t_c * pow(r_j, -0.75); if (Flags.debug[GEO]) { printf("Ring : %ld \n", i); printf("alpha : %f \n", alpha); printf("height : %f \n", ring_height); printf("temp : %f \n", t_temp); } /* ----------------- loop over phi -------------------------------- */ for (j = 0; j < n_phi; ++j) { phi = ((2.0*PI/n_phi)*j)-Phase; /* MPR 23.05.2002 */ phi += dphi; sin_phi = sin(phi); cos_phi = cos(phi); /* PR added 07.05.02 is needed for eclipsing * by disk verification */ SurfPtr->RadiusIn = r_i_old; SurfPtr->RadiusOut = r_i_new; /* phi counted from x to y axis */ SurfPtr->lambda = (float) r_j * cos_phi; SurfPtr->mu = (float) r_j * sin_phi; /* z is thickness/2.0 + (H/R * r_j)/2.0 */ SurfPtr->nu = (float) ring_height; /* for a flaring disk, the normal vector points up and inwards, * so we need to switch sign for l, m */ SurfPtr->l = (float) (-cos_phi) * sin_alpha; if (fabs(SurfPtr->l) <= DBL_EPSILON) { SurfPtr->l = 0.0; } SurfPtr->m = (float) (-sin_phi) * sin_alpha; if (fabs(SurfPtr->m) <= DBL_EPSILON) { SurfPtr->m = 0.0; } SurfPtr->n = (float) cos_alpha; if (fabs(SurfPtr->n) <= DBL_EPSILON) { SurfPtr->n = 0.0; } #if 0 if (Flags.warp == ON) { TileDisk1( warpfactor[i], SurfPtr ); TileDisk2( warpfactor[i], SurfPtr ); } if (tilt_flag) { TileDisk1( TiltAngle, SurfPtr ); TileDisk2( TiltAngle, SurfPtr ); } /* Calculate the layer vectors */ CalculateLayerVecE1( SurfPtr ); CalculateLayerVecE2( SurfPtr ); #endif SurfPtr->AreaType = TOP_SEGMENT; /* Angle between x-Axis and Middle of the Element */ SurfPtr->AreaAngle = j*360/n_phi; SurfPtr->MyNum = NumElem; SurfPtr->rho = 1.00; /* ???*/ SurfPtr->SpotNum = 0; SurfPtr->rad = r_j; /* Is needed for eclipse test */ SurfPtr->area = DSurfCor; Surf = Surf + DSurfCor; /* --------------- INITIALIZE SURFACE TEMP ----------------- */ SurfPtr->temp = t_temp; hs_phi = fmod(phi, 2.0 * PI); hs_phi = (hs_phi < 0) ? (2.0 * PI + hs_phi) : hs_phi; if ( (DIST_CIRCLE(hs_phi, hs_long) <= hs_ext) && (r_j >= (Rout - hs_depth) )) { SurfPtr->temp = hs_temp; SurfPtr->SpotNum = 1; } ++NumElem; /* increment surface element count */ ++SurfPtr; /* some elements of the structure are not filled at the moment */ } /* end of loop over phi */ r_i_old = r_i_new; /* end of one ring */ } /* end of loop over rings (n_r) */ /**********************************************************************/ /* */ /* ---------------------- Calculate outer side of disk ------------- */ /* */ /**********************************************************************/ /* area is outer side divided by rings/segments */ DSurf = ((2.0*PI*Rout*Thick_out) / n_phi) / BinPtr->RingsOut; /* vertical offset unit (half ring height) */ offset_out = (Thick_out/BinPtr->RingsOut)/2.0; /* ----------------- loop over rings (n_r) ---------------------- */ for (i = 0; i < BinPtr->RingsOut ; ++i) { phi = 0.0; /* Correct ?????*/ /* Height of this ring. Start at top, * subtract half ring height, subtract (ring height * #ring) */ ring_height = ( (diskAdd + diskProp * pow(Rout, diskExp)) / 2.0 ) - offset_out - i * (Thick_out/BinPtr->RingsOut); if (Flags.debug[GEO]) { printf("OuterRing %03ld: %f\n", i, ring_height); } /* check if enough memory */ if ((NumElem+n_phi) >= (int) MaximumElements) { if (Flags.debug[GEO]) printf("Outer side of disk: %ld >= %d\n", (NumElem+n_phi), (int) MaximumElements); #ifdef _WITH_GTK if (Flags.interactive == OFF) nf_error(_(errmsg[10])); else { if (Flags.WantFit == OFF && Flags.WantMap == OFF) make_dialog(_(errmsg[10])); return(8); } #else nf_error(_(errmsg[10])); #endif } /* dphi = i * ((2.0*PI/n_phi)/BinPtr->RingsOut); */ /* shift rings in phi by random amount to reduce numerical artifacts */ #ifdef _WITH_OPENGL if (Flags.use_opengl == OFF) #endif dphi = (2.0*PI * (rand() / (RAND_MAX + 1.0))); if (t_flag < 2) /* isothermal */ t_temp = T1; else t_temp = t_c * pow(Rout, -0.75); /* ----------------- loop over phi -------------------------------- */ for (j = 0; j < n_phi; ++j) { phi = ((2.0*PI/n_phi)*j)-Phase; /* MPR 23.05.2002 */ /* Shift rings in phi to avoid spikes in light curve. */ phi += dphi; sin_phi = sin(phi); cos_phi = cos(phi); SurfPtr->RadiusIn = Rout; /* PR added 07.05.02 is needed for */ SurfPtr->RadiusOut = Rout; /* eclipsing by disk verification */ SurfPtr->lambda = (float) Rout * cos_phi; SurfPtr->mu = (float) Rout * sin_phi; SurfPtr->nu = (float) ring_height; SurfPtr->l = (float) cos_phi; if (fabs(SurfPtr->l) <= DBL_EPSILON){ SurfPtr->l = 0.0; } SurfPtr->m = (float) sin_phi; // PR -sin -> sin 21.05.02 if (fabs(SurfPtr->m) <= DBL_EPSILON){ SurfPtr->m = 0.0; } SurfPtr->n = 0.0; #if 0 if (tilt_flag) { TileDisk1(TiltAngle,SurfPtr); TileDisk2(TiltAngle,SurfPtr); } /* Calculate the layer vectors */ CalculateLayerVecE1(SurfPtr); CalculateLayerVecE2(SurfPtr); #endif SurfPtr->AreaType = OUT_RECTANGLE; /* Angle between x-Axis and Middle of the Element */ SurfPtr->AreaAngle = (j*360.0)/n_phi; SurfPtr->AreaHeight = Thick_out/BinPtr->RingsOut; SurfPtr->AreaWidth = (2*PI*Rout)/n_phi; SurfPtr->MyNum = NumElem; SurfPtr->ring = i; SurfPtr->rad = Rout; /* Is needed for eclipse test */ SurfPtr->area = DSurf; Surf = Surf + DSurf; /* --------------- INITIALIZE SURFACE TEMP ----------------- */ SurfPtr->temp = t_temp; hs_phi = fmod(phi, 2.0 * PI); hs_phi = (hs_phi < 0) ? (2.0 * PI + hs_phi) : hs_phi; if ( (DIST_CIRCLE(hs_phi, hs_long) <= hs_ext) ) { SurfPtr->temp = hs_temp; SurfPtr->SpotNum = 1; } ++NumElem; /* increment surface element count */ ++SurfPtr; /* some elements of the structure are not filled at the moment */ } /* end of loop over phi */ } /* end of loop over rings (n_r) */ /**********************************************************************/ /* */ /* ---------------------- Calculate bottom of disk ----------------- */ /* */ /**********************************************************************/ /* Here we loop from outer to inner border */ r_outer_sqr = (r_i_old*r_i_old); /* Area of an element. This can be lifted outside the loop, since * the progression of radii ensures equal areas for each element. */ DSurf = (PI*P_0)/n_phi; /* ----------------- loop over rings (n_r) ---------------------- */ for (i = 0; i < n_r; ++i) { phi = 0.0; /* check if enough memory */ if ((NumElem+n_phi) >= (int) MaximumElements) { if (Flags.debug[GEO]) printf("Bottom of disk: %ld >= %d\n", (NumElem+n_phi), (int) MaximumElements); #ifdef _WITH_GTK if (Flags.interactive == OFF) nf_error(_(errmsg[10])); else { if (Flags.WantFit == OFF && Flags.WantMap == OFF) make_dialog(_(errmsg[10])); return(8); } #else nf_error(_(errmsg[10])); #endif } r_i_new = sqrt(r_outer_sqr - (i+1)*P_0); /* radius of next circle */ r_j = (r_i_new+r_i_old)/2.0; /* Center of the element */ warp_no = (int) n_r-i; ring_height = -1.0 * ((diskAdd + diskProp * pow(r_j, diskExp)) / 2.0); /* angle of rise for this ring */ h_i_old = -1.0 * ( (diskAdd + diskProp * pow(r_i_old, diskExp)) / 2.0); h_i_new = -1.0 * ( (diskAdd + diskProp * pow(r_i_new, diskExp)) / 2.0); rise_angle = atan ( fabs(h_i_old - h_i_new) / fabs(r_i_old - r_i_new) ); /* Corrected surface area */ DSurfCor = DSurf / cos ( rise_angle ); /* shift rings in phi by random amount to reduce numerical artifacts */ #ifdef _WITH_OPENGL if (Flags.use_opengl == OFF) #endif dphi = (2.0*PI * (rand() / (RAND_MAX + 1.0))); /* Angle of rise from the middle to the edge of the Accretion Disk */ if (t_flag == 0) alpha = (- 0.5 * diskProp); else alpha = -rise_angle /* (- 0.5 * diskProp) * diskExp * pow(r_j, (diskExp - 1.0)) */; sin_alpha = sin(alpha); cos_alpha = cos(alpha); if (t_flag < 2) /* isothermal */ t_temp = T1; else t_temp = t_c * pow(r_j, -0.75); /* ----------------- loop over phi -------------------------------- */ for (j = 0; j < n_phi; ++j) { phi=((2.0*PI/n_phi)*j)-Phase; /* MPR 23.05.2002 */ phi += dphi; sin_phi = sin(phi); cos_phi = cos(phi); /* PR added 07.05.02 is needed for eclipsing * by disk verification * BEWARE!!! r_i_new und r_i_old muessen hier vertauscht sein * da die Schleife jetzt von Aussen nach Innen laeuft!!!! */ SurfPtr->RadiusIn = r_i_new; SurfPtr->RadiusOut = r_i_old; SurfPtr->lambda = (float) r_j * cos_phi; SurfPtr->mu = (float) r_j * sin_phi; /* z negative (this is lower surface) */ SurfPtr->nu = (float) ring_height; SurfPtr->l = (float) (-cos_phi) * sin_alpha; if (fabs(SurfPtr->l) <= DBL_EPSILON) { SurfPtr->l = 0.0; } SurfPtr->m = (float) (-sin_phi) * sin_alpha; if (fabs(SurfPtr->m) <= DBL_EPSILON) { SurfPtr->m = 0.0; } SurfPtr->n = (float) -1.0 * (cos_alpha); if (fabs(SurfPtr->n) <= DBL_EPSILON) { SurfPtr->n = 0.0; } #if 0 if (Flags.warp == ON) { TileDisk1(warpfactor[warp_no],SurfPtr); TileDisk2(warpfactor[warp_no],SurfPtr); } if (tilt_flag) { TileDisk1(TiltAngle,SurfPtr); TileDisk2(TiltAngle,SurfPtr); } /* Calculate the layer vectors */ CalculateLayerVecE1(SurfPtr); CalculateLayerVecE2(SurfPtr); #endif SurfPtr->AreaType = BOTTOM_SEGMENT; /* Angle between x-Axis and Middle of the Element */ SurfPtr->AreaAngle = j*360/n_phi; SurfPtr->MyNum = NumElem; SurfPtr->ring = (n_r - i) - 1; SurfPtr->rad = r_j; /* Is needed for eclipse test */ SurfPtr->area = DSurfCor; Surf = Surf + DSurfCor; /* --------------- INITIALIZE SURFACE TEMP ----------------- */ SurfPtr->temp = t_temp; hs_phi = fmod(phi, 2.0 * PI); hs_phi = (hs_phi < 0) ? (2.0 * PI + hs_phi) : hs_phi; if ( (DIST_CIRCLE(hs_phi, hs_long) <= hs_ext) && (r_j >= (Rout - hs_depth) )) { SurfPtr->temp = hs_temp; SurfPtr->SpotNum = 1; } ++NumElem; /* increment surface element count */ ++SurfPtr; /* some elements of the structure are not filled at the moment */ } /* end of loop over phi */ r_i_old = r_i_new; /* end of one ring */ } /* end of loop over rings (n_r) */ /**********************************************************************/ /* */ /* ---------------------- Calculate inner side of disk ------------- */ /* */ /**********************************************************************/ /* area is inner side divided by rings/segments */ DSurf = ((2.0*PI*Rin*Thick_in)/n_phi) / BinPtr->RingsIn; /* vertical offset unit (half ring height) */ offset_in = (Thick_in/BinPtr->RingsIn)/2.0; /* ----------------- loop over rings (n_r) ---------------------- */ for (i = 0; i < BinPtr->RingsIn ; ++i) { phi = 0.0; /* Correct ?????*/ /* Height of this ring. Start at bottom, * add half ring height, add (ring height * #ring) */ ring_height = -1.0 * ( (diskAdd + diskProp * pow(Rin, diskExp)) / 2.0 ) + offset_in + i * (Thick_in/BinPtr->RingsIn); if (Flags.debug[GEO]) { printf("InnerRing %0ld: %f\n", i, ring_height); } /* check if enough memory */ if ((NumElem+n_phi) >= (int) MaximumElements) { if (Flags.debug[GEO]) printf("Inner side of disk: %ld >= %d\n", (NumElem+n_phi), (int) MaximumElements); #ifdef _WITH_GTK if (Flags.interactive == OFF) nf_error(_(errmsg[10])); else { if (Flags.WantFit == OFF && Flags.WantMap == OFF) make_dialog(_(errmsg[10])); return(8); } #else nf_error(_(errmsg[10])); #endif } /* dphi = i * ((2.0*PI/n_phi)/BinPtr->RingsIn); */ /* shift rings in phi by random amount to reduce numerical artifacts */ #ifdef _WITH_OPENGL if (Flags.use_opengl == OFF) #endif dphi = (2.0*PI * (rand() / (RAND_MAX + 1.0))); if (t_flag < 2) /* isothermal */ t_temp = T1; else t_temp = t_c * pow(Rin, -0.75); /* ----------------- loop over phi -------------------------------- */ for (j = 0; j < n_phi; ++j) { phi = ((2.0*PI/n_phi)*j)-Phase; /* MPR 23.05.2002 */ /* Shift rings in phi to avoid spikes in light curve. */ phi += dphi; sin_phi = sin(phi); cos_phi = cos(phi); SurfPtr->RadiusIn = Rin; /* PR added 07.05.02 is needed for */ SurfPtr->RadiusOut = Rin; /* eclipsing by disk verification */ SurfPtr->lambda = (float) Rin * cos_phi; SurfPtr->mu = (float) Rin * sin_phi; SurfPtr->nu = (float) ring_height; /* Invert sign of cos(phi),/sin(phi), since the normal vector * points towards the z-axis here. */ SurfPtr->l = (float) -1.0 * (cos_phi); if (fabs(SurfPtr->l) <= DBL_EPSILON) { SurfPtr->l = 0.0; } SurfPtr->m = (float) -1.0 * (sin_phi); if (fabs(SurfPtr->m) <= DBL_EPSILON) { SurfPtr->m = 0.0; } SurfPtr->n = 0.0; #if 0 if (tilt_flag) { TileDisk1(TiltAngle,SurfPtr); TileDisk2(TiltAngle,SurfPtr); } /* Calculate the layer vectors */ CalculateLayerVecE1(SurfPtr); CalculateLayerVecE2(SurfPtr); #endif SurfPtr->AreaType = IN_RECTANGLE; /* Angle between x-Axis and Middle of the Element */ SurfPtr->AreaAngle = j*360/n_phi; SurfPtr->AreaHeight = Thick_in/BinPtr->RingsIn; SurfPtr->AreaWidth = (2*PI*Rin)/n_phi; SurfPtr->MyNum = NumElem; SurfPtr->ring = i; SurfPtr->rad = Rin; /* Is needed for eclipse test */ SurfPtr->area = DSurf; Surf = Surf + DSurf; /* --------------- INITIALIZE SURFACE TEMP ----------------- */ SurfPtr->temp = t_temp; hs_phi = fmod(phi, 2.0 * PI); hs_phi = (hs_phi < 0) ? (2.0 * PI + hs_phi) : hs_phi; if ( (DIST_CIRCLE(hs_phi, hs_long) <= hs_ext) && (Rin >= (Rout - hs_depth) )) { SurfPtr->temp = hs_temp; SurfPtr->SpotNum = 1; } ++NumElem; /* increment surface element count */ ++SurfPtr; /* some elements of the structure are not filled at the moment */ } /* end of loop over phi */ } /* end of loop over rings (n_r) */ BinPtr->NumElem = NumElem; /* Number of elements*/ if (Flags.debug[VERBOSE] == ON) { printf("\n"); printf(_(" Divide Disk %3d : %6ld Surface elements\n"), Comp+1, BinPtr->NumElem); } BinPtr->Surface = Surf; /*Added 29-04-2002, EG/MPR */ /* ---------------------------- mean temperature ------------------- */ Weight = 0.0; TempMean = 0.0; SurfPtr = Surface[Comp]; for (i = 0; i < BinPtr->NumElem; ++i) { Area = SurfPtr->area; TempMean = TempMean + SurfPtr->temp * Area; Weight = Weight + Area; /*------- Sollte hier nicht irgendwo ein SurfPtr++ stehen????? */ /* eigentlich schon ... RW Tue Sep 11 12:54:20 CEST 2001 */ ++SurfPtr; } BinPtr->TempMean = TempMean/Weight; /**********************************************************************/ /* */ /* Flux intensity computation */ /* */ /**********************************************************************/ if (Flags.blackbody == ON) { BlackbodyFlux(Disk); } else { Binary[Disk].log_g = 3.5; ModelFlux(Disk); } /* ---------------- output --------------------- */ if (Flags.debug[VERBOSE] == ON) { /* printf("\n Star: %5d\n", (Comp+1)); */ printf(_(" Surface Area, Mean Gravity (dimensionless): %8.6f %8.6f\n"), BinPtr->Surface, BinPtr->Gravity); printf(_(" Temperature is %8.6f\n"), BinPtr->Temperature); printf("\n"); } return(0); } #endif