/* 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 <stdlib.h>
#include <string.h>
#include <float.h>
#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
syntax highlighted by Code2HTML, v. 0.9.1