/* 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"
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Simple reflection treatment (point source)
@tip Includes penumbral correction
@param (int) Comp The stellar component
@return (void)
@heading Light Curve
****************************************************************************/
void SimpleReflect(int Comp)
{
double penumbral_corr, penumbral_h; /* penumbral correction */
double penumbral_a, penumbral_alpha; /* penumbral correction */
double penumbral_sin; /* penumbral correction */
double cosalpha, sinalpha, sumsin; /* reflection angle */
double lvec, mvec, nvec; /* surface normal */
double T1, T2, R2, Tr, rho; /* blackbody */
long NumElem; /* loop counter */
SurfaceElement *SurfPtr; /* surface pointer */
BinaryComponent *BinPtr; /* star pointer */
T1 = Binary[Comp+1].Temperature;
SurfPtr = Surface[Comp];
BinPtr = &Binary[Comp];
if (Binary[Primary].RocheFill >= 0.3 || Binary[Secondary].RocheFill >= 0.3)
WARNING(_("Simple Reflection not recommended for large fill factors"));
if (Comp == 0) {
T2 = Binary[Comp+1].Temperature;
R2 = Binary[Comp+1].Radius;
} else {
T2 = Binary[Comp-1].Temperature;
R2 = Binary[Comp-1].Radius;
}
for (NumElem = 0; NumElem < BinPtr->NumElem; ++NumElem) {
rho = SurfPtr->rho;
lvec = SurfPtr->l;
mvec = SurfPtr->m;
nvec = SurfPtr->n;
/* Eq. 1-19, Djurasevic (1992). Lambda, mu, nu are already
* multiplied with radius.
*/
cosalpha = rho *
( lvec * ( 1.0-SurfPtr->lambda)
- SurfPtr->mu * mvec
- SurfPtr->nu * nvec);
sinalpha = sqrt(1.0 - SQR(cosalpha));
/* Following Eq. 1-19, Djurasevic (1992),
* rho* := 1/rho is distance to neighbouring star centre.
* Radius * rho = radius / rho* is approximately sin of angle
* subtended by radius of neighbouring star.
*/
penumbral_sin = R2*rho;
penumbral_sin = (penumbral_sin > 1.0) ? 1.0 : penumbral_sin;
penumbral_corr = sqrt(1.0-SQR(penumbral_sin));
/* This is the sinus of the angle subtended by the visible sector of the
* neighbouring star (vertical to the horizon). Approximately.
*/
sumsin = cosalpha + penumbral_sin;
if (cosalpha <= (-DBL_EPSILON) && sumsin <= (-DBL_EPSILON)) {
/* fully invisible, do nothing */
;
#if 0
if (Comp == 0) {
printf("NADA %05ld %12.6f %12.6f %12.6f %12.6f\n",
NumElem, SurfPtr->eta, cosalpha, penumbral_sin, sumsin);
}
#endif
}
else if (cosalpha >= 0 && sumsin >= (2.0 * penumbral_sin)) {
SurfPtr->AreaType = FRONTSIDE;
/* fully visible */
T1 = SurfPtr->temp;
Tr = T2/T1; Tr = Tr * Tr * Tr * Tr;
SurfPtr->temp = T1 *
sqrt(sqrt(1.0 + BinPtr->Albedo * 0.5 * cosalpha *
(1.0 - sqrt(1.0 - SQR(penumbral_sin)))*Tr));
#if 0
if (Comp == 0) {
printf("FULL %05ld %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f\n",
NumElem, SurfPtr->eta, cosalpha, penumbral_sin, sumsin,
T1, Tr, SurfPtr->temp);
}
#endif
}
else {
SurfPtr->AreaType = FRONTSIDE;
T1 = SurfPtr->temp;
Tr = T2/T1; Tr = Tr * Tr * Tr * Tr;
/* partially visible, apply penumbral correction */
/* The height of the visible disc sector in units of the radius;
* sumsin is sin(visible height), penumbral_sin is sin(radius).
*/
penumbral_h = sumsin/penumbral_sin;
if (cosalpha >= 0.0) {
/* We see more than half the star, so we compute what we need to
* subtract from the full disk area (2 PI). We are interested in
* the area of the invisible segment of height (2.0 - penumbral_h).
*/
penumbral_corr = PI + PI; /* set to full disk */
penumbral_h = 2.0 - penumbral_h;
/* we calculate the segment of a disk, r=1.0 */
/* penumbral_a is **half** the chord of the segment */
penumbral_a = sqrt(penumbral_h+penumbral_h - SQR(penumbral_h));
/* Okt 14 15:38:27 MEST 1999 -- fix rounding error */
penumbral_a = MIN(1.0, penumbral_a);
/* penumbral_alpha is half the arc length of the segment */
penumbral_alpha = asin( CLAMP(penumbral_a, 0.0, 1.0) );
/* segment size is:
* 0.5*(arc_length - chord_length * (1.0 - segment_height))
* (penumbral_alpha - penumbral_a * (1.0 - penumbral_h ))
*/
/* subtract segment */
penumbral_corr = penumbral_corr
- (penumbral_alpha - penumbral_a * (1.0 - penumbral_h));
/* Convert to fractional area for correction factor */
penumbral_corr = penumbral_corr / (2.0 * PI);
/* Eq. 1-23, Djurasevic (1992). Here, penumbral_corr is the
* computed fractional visible disk. Cosalpha is increased
* by 0.5 * angle subtended by neighbouring star vertical
* to the horizon.
*/
SurfPtr->temp = T1 *
sqrt(sqrt(1.0 +
penumbral_corr * BinPtr->Albedo * 0.5
* (cosalpha + (0.5 * sumsin)) *
(1.0 - sqrt(1.0 - SQR(penumbral_sin)))*Tr));
} else {
/* We see less than half the star, so the visible part is
* a segment of a circle.
*/
/* penumbral_a is **half** the chord of the segment */
penumbral_a = sqrt(2.0*penumbral_h - SQR(penumbral_h));
/* Okt 14 15:38:27 MEST 1999 -- fix rounding error */
penumbral_a = MIN(1.0, penumbral_a);
/* penumbral_alpha is half the arc length of the segment */
penumbral_alpha = asin( CLAMP(penumbral_a, 0.0, 1.0) );
/* segment size is:
* 0.5*(arc_length - chord_length * (1.0 - segment_height))
* (penumbral_alpha - penumbral_a * (1.0 - penumbral_h ))
*/
penumbral_corr =
penumbral_alpha - penumbral_a * (1.0 - penumbral_h);
/* Convert to fractional area for correction factor */
penumbral_corr = penumbral_corr / (2.0 * PI);
/* Eq. 1-23, Djurasevic (1992). Here, penumbral_corr is the
* computed fractional visible disk. Cosalpha is replaced
* by 0.5 * angle subtended by neighbouring star vertical
* to the horizon.
*/
SurfPtr->temp = T1 *
sqrt(sqrt(1.0 +
penumbral_corr * BinPtr->Albedo * 0.5 *
(0.5 * sumsin) *
(1.0 - sqrt(1.0 - SQR(penumbral_sin)))*Tr));
}
#if 0
if (Comp == 0) {
printf("PART %05ld %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f %12.6f\n",
NumElem, SurfPtr->eta, cosalpha, penumbral_sin, sumsin,
T1, Tr, SurfPtr->temp);
}
#endif
}
if (SurfPtr->temp > DBL_MAX || SurfPtr->temp < DBL_MIN)
SurfPtr->temp = T1;
++SurfPtr;
}
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Detailed reflection treatment (mutual irradiation)
@param (void)
@return (void)
@heading Light Curve
****************************************************************************/
void LightReflect()
{
#define NSQR_TAB 1000
static double SumDelta_[2][MAXELEMENTS]; /* delta flux */
static double Temp_[2][MAXELEMENTS]; /* temperature */
static double sqtab[NSQR_TAB+1] = { 0.0 };/* sqrt lookup table */
static double u1[2], u2[2]; /* limb darkening */
static double NormLimb_[2]; /* normalization */
float *lu1, *lu2, *lT; /* pointer to LDC table */
long i, j, k; /* loop counter */
long NumLarge, NumSmall; /* # of surface elements */
int LargeStar, SmallStar; /* the stars */
double *Init_ptr; /* used to init SumDelta_ */
double dx, dy, dz; /* distance */
double cosgamma_ji, cosgamma_ij; /* viewing angles */
double distance_ij; /* distance */
double Delta_F, Delta_Ii, Delta_Ij; /* deltas */
double Tn, Tsmall, Tlarge,Ttmp; /* temperatures */
int Nmax, Nnow; /* where do we stand ? */
SurfaceElement *SurfPtrS; /* surface pointer */
SurfaceElement *SurfPtrL; /* surface pointer */
int LoopHigh, LoopLow; /* upper/lower limit for loop */
#ifdef _WITH_GTK
char message[64]; /* statusbar message */
#endif
int gsize = numprocs; /* # processes */
#if defined(_WITH_MPI_COARSE)
static double RecvBuf[MAXELEMENTS]; /* delta flux */
int *displs = NULL; /* displacements */
int *rcount = NULL; /* # received data */
int *sdispls = NULL; /* displacements */
int *srcount = NULL; /* # received data */
int stride; /* # data/process */
#endif
/* Bolometric limb darkening coefficients log_g=4.0 sqrt law */
/* Claret A. 2000, Astron. Astrophys. 363, 1081 */
/* PHOENIX */
static
float TH[35] ={ 2000., 2200., 2400., 2600., 2800., 3000., 3200., 3400., 3600., 3800.,
4000., 4200., 4400., 4600., 4800., 5000., 5200., 5400., 5600., 5800.,
7000., 7200., 7400., 7600., 7800., 8000., 8200., 8400., 8600., 8800.,
9000., 9200., 9400., 9600., 9800.};
static
float hu1_35[35] ={0.3287, 0.3287, 0.1738, -.0272, -.1646, -.2050, -.0797, 0.2362, 0.4438, 0.4517,
0.0949, 0.0834, 0.0825, 0.0882, 0.0891, 0.0736, 0.0627, 0.0504, 0.0409, 0.0347,
0.0070, -.0008, -.0021, -.0110, -.0202, -.0293, -.0373, -.0467, -.0520, -.0575,
-.0659, -.0696, -.0894, -.0864, -.0933};
static
float hu2_35[35] ={0.4955, 0.5954, 0.7281, 0.8824, 0.9538, 0.9447, 0.7892, 0.4551, 0.2474, 0.2747,
0.5918, 0.5662, 0.5492, 0.5209, 0.5043, 0.5188, 0.5141, 0.5102, 0.5032, 0.4928,
0.4439, 0.4433, 0.4360, 0.4363, 0.4389, 0.4423, 0.4451, 0.4523, 0.4514, 0.4509,
0.4640, 0.4612, 0.4871, 0.4763, 0.4844};
static
float hu1_40[35] ={0.4206, 0.3032, 0.1503, 0.0008, -.1148, -.2062, -.1686, 0.0298, 0.2389, 0.3105,
0.1565, 0.0862, 0.0423, 0.0319, 0.0244, 0.0337, 0.0286, 0.0169, 0.0053, -.0028,
-.0230, -.0360, -.0360, -.0437, -.0522, -.0604, -.0639, -.0661, -.0703, -.0780,
-.0831, -.0892, -.0888, -.0867, -.0882};
static
float hu2_40[35] ={0.5268, 0.6257, 0.7510, 0.8575, 0.9050, 0.9289, 0.8508, 0.6266, 0.4275, 0.3862,
0.5572, 0.6300, 0.6649, 0.6582, 0.6528, 0.5947, 0.5830, 0.5791, 0.5759, 0.5684,
0.5132, 0.5186, 0.5093, 0.5075, 0.5080, 0.5090, 0.5050, 0.5067, 0.4985, 0.4951,
0.4907, 0.4978, 0.4919, 0.4810, 0.4787};
static
float hu1_45[35] ={0.3977, 0.3054, 0.1673, 0.0103, -.0987, -.1867, -.1923, -.1072, 0.0039, 0.0986,
0.1549, 0.1639, 0.0658, 0.0498, 0.0442, 0.0487, 0.0467, 0.0387, 0.0285, 0.0202,
-.0062, -.0119, -.0183, -.0252, -.0254, -.0328, -.0404, -.0473, -.0506, -.0534,
-.0614, -.0639, -.0652, -.0651, -.0679};
static
float hu2_45[35] ={0.5548, 0.6343, 0.7467, 0.8570, 0.8996, 0.9134, 0.8701, 0.7447, 0.6249, 0.5583,
0.5273, 0.5568, 0.6114, 0.6045, 0.5934, 0.5625, 0.5494, 0.5432, 0.5394, 0.5332,
0.4856, 0.4820, 0.4795, 0.4782, 0.4701, 0.4681, 0.4678, 0.4670, 0.4628, 0.4627,
0.4535, 0.4493, 0.4439, 0.4371, 0.4370};
static
float hu1_50[35] ={0.3768, 0.2876, 0.1685, 0.0287, -.0794, -.1494, -.1966, -.1713, -.1069, -.0790,
-.0334, 0.0419, 0.0870, 0.0555, 0.0501, 0.0480, 0.0466, 0.0427, 0.0312, 0.0236,
-.0007, -.0053, -.0105, -.0160, -.0223, -.0243, -.0298, -.0368, -.0438, -.0493,
-.0526, -.0561, -.0582, -.0549, -.0610};
static
float hu2_50[35] ={0.5694, 0.6540, 0.7555, 0.8574, 0.9003, 0.8997, 0.8844, 0.8121, 0.7133, 0.6935,
0.6673, 0.6255, 0.5991, 0.5971, 0.5764, 0.5589, 0.5461, 0.5366, 0.5362, 0.5303,
0.4794, 0.4743, 0.4701, 0.4669, 0.4648, 0.4586, 0.4557, 0.4546, 0.4541, 0.4516,
0.4439, 0.4400, 0.4325, 0.4211, 0.4226};
/* ATLAS */
static
float TK[76] ={ 3500., 3750., 4000., 4250., 4500., 4750., 5000., 5250., 5500.,
5750., 6000., 6250., 6500., 6750., 7000., 7250., 7500., 7750.,
8000., 8250., 8500., 8750., 9000., 9250., 9500., 9750., 10000.,
10250., 10500., 10750., 11000., 11250., 11500., 11750., 12000., 12250.,
12500., 12750., 13000., 14000., 15000., 16000., 17000., 18000., 19000.,
20000., 21000., 22000., 23000., 24000., 25000., 26000., 27000., 28000.,
29000., 30000., 31000., 32000., 33000., 34000., 35000., 36000., 37000.,
38000., 39000., 40000., 41000., 42000., 43000., 44000., 45000., 46000.,
47000., 48000., 49000., 50000.};
/* data for 3.5 only until 31000K, repeat last data point to 50000K */
static
float lu1_35[76] ={0.0195, 0.1686, 0.2396, 0.2580, 0.2647, 0.2517, 0.2214, 0.1848, 0.1494,
0.1187, 0.0942, 0.0787, 0.0695, 0.0637, 0.0651, 0.0909, 0.1712, 0.2590,
0.3497, 0.3857, 0.4060, 0.4419, 0.4702, 0.4893, 0.4971, 0.4983, 0.4939,
0.4884, 0.4829, 0.4789, 0.4778, 0.4798, 0.4865, 0.4936, 0.5025, 0.5116,
0.5211, 0.5282, 0.5355, 0.5585, 0.5684, 0.5741, 0.5770, 0.5770, 0.5748,
0.5666, 0.5522, 0.5341, 0.5096, 0.4796, 0.4445, 0.4058, 0.3654, 0.3332,
0.3172, 0.3163, 0.3221, 0.3221, 0.3221, 0.3221, 0.3221, 0.3221, 0.3221,
0.3221, 0.3221, 0.3221, 0.3221, 0.3221, 0.3221, 0.3221, 0.3221, 0.3221,
0.3221, 0.3221, 0.3221, 0.3221};
static
float lu2_35[76] ={0.6125, 0.4919, 0.4355, 0.4252, 0.4212, 0.4393, 0.4754, 0.5153, 0.5510,
0.5795, 0.5997, 0.6119, 0.6207, 0.6268, 0.6273, 0.5970, 0.5159, 0.4318,
0.3332, 0.2941, 0.2743, 0.2366, 0.2119, 0.2004, 0.2027, 0.2109, 0.2239,
0.2367, 0.2481, 0.2564, 0.2601, 0.2594, 0.2523, 0.2450, 0.2353, 0.2258,
0.2159, 0.2095, 0.2029, 0.1842, 0.1803, 0.1788, 0.1781, 0.1792, 0.1812,
0.1896, 0.2045, 0.2225, 0.2473, 0.2777, 0.3133, 0.3513, 0.3890, 0.4141,
0.4202, 0.4118, 0.4003, 0.4003, 0.4003, 0.4003, 0.4003, 0.4003, 0.4003,
0.4003, 0.4003, 0.4003, 0.4003, 0.4003, 0.4003, 0.4003, 0.4003, 0.4003,
0.4003, 0.4003, 0.4003, 0.4003};
/* data for 4.0 only until 39000K, repeat last data point to 50000K */
static
float lu1_40[76] ={-.0826, 0.0527, 0.1948, 0.2525, 0.2628, 0.2528, 0.2270, 0.1958, 0.1612,
0.1284, 0.1048, 0.0894, 0.0823, 0.0758, 0.0703, 0.0687, 0.0773, 0.1440,
0.2357, 0.3501, 0.4075, 0.4410, 0.4783, 0.5028, 0.5159, 0.5194, 0.5158,
0.5078, 0.4976, 0.4899, 0.4854, 0.4830, 0.4861, 0.4910, 0.4979, 0.5065,
0.5153, 0.5239, 0.5312, 0.5577, 0.5714, 0.5758, 0.5765, 0.5754, 0.5723,
0.5643, 0.5517, 0.5356, 0.5109, 0.4817, 0.4467, 0.4067, 0.3633, 0.3201,
0.2833, 0.2593, 0.2509, 0.2510, 0.2555, 0.2669, 0.2895, 0.3283, 0.3643,
0.3895, 0.4050, 0.4050, 0.4050, 0.4050, 0.4050, 0.4050, 0.4050, 0.4050,
0.4050, 0.4050, 0.4050, 0.4050};
static
float lu2_40[76] ={0.6764, 0.5939, 0.4773, 0.4296, 0.4221, 0.4372, 0.4691, 0.5037, 0.5394,
0.5709, 0.5919, 0.6044, 0.6085, 0.6158, 0.6233, 0.6252, 0.6154, 0.5483,
0.4588, 0.3313, 0.2653, 0.2295, 0.1926, 0.1754, 0.1730, 0.1814, 0.1961,
0.2141, 0.2333, 0.2478, 0.2571, 0.2629, 0.2613, 0.2574, 0.2507, 0.2420,
0.2333, 0.2250, 0.2187, 0.1970, 0.1911, 0.1946, 0.2000, 0.2048, 0.2102,
0.2202, 0.2345, 0.2519, 0.2787, 0.3095, 0.3460, 0.3871, 0.4303, 0.4704,
0.4998, 0.5123, 0.5076, 0.4955, 0.4825, 0.4626, 0.4302, 0.3806, 0.3393,
0.3144, 0.3012, 0.3012, 0.3012, 0.3012, 0.3012, 0.3012, 0.3012, 0.3012,
0.3012, 0.3012, 0.3012, 0.3012};
/* data for 4.5 only until 49000K, repeat last data point to 50000K */
static
float lu1_45[76] ={-.1460, -.0807, 0.0616, 0.2086, 0.2562, 0.2501, 0.2286, 0.1974, 0.1653,
0.1376, 0.1133, 0.0994, 0.0941, 0.0887, 0.0833, 0.0804, 0.0788, 0.0877,
0.1198, 0.2093, 0.3326, 0.4272, 0.4742, 0.5094, 0.5321, 0.5401, 0.5382,
0.5291, 0.5186, 0.5087, 0.4970, 0.4912, 0.4911, 0.4933, 0.4980, 0.5046,
0.5116, 0.5189, 0.5261, 0.5513, 0.5676, 0.5749, 0.5753, 0.5715, 0.5647,
0.5567, 0.5431, 0.5249, 0.5002, 0.4699, 0.4320, 0.3916, 0.3458, 0.2990,
0.2551, 0.2184, 0.1970, 0.1880, 0.1887, 0.1962, 0.2055, 0.2247, 0.2591,
0.3008, 0.3380, 0.3648, 0.3824, 0.3944, 0.4054, 0.4158, 0.4268, 0.4366,
0.4457, 0.4537, 0.4611, 0.4611};
static
float lu2_45[76] ={0.7201, 0.6991, 0.5982, 0.4710, 0.4286, 0.4394, 0.4668, 0.5024, 0.5357,
0.5623, 0.5841, 0.5948, 0.5965, 0.6024, 0.6097, 0.6144, 0.6195, 0.6136,
0.5828, 0.4918, 0.3527, 0.2412, 0.1905, 0.1584, 0.1443, 0.1489, 0.1644,
0.1864, 0.2080, 0.2269, 0.2467, 0.2580, 0.2611, 0.2612, 0.2579, 0.2522,
0.2461, 0.2397, 0.2334, 0.2136, 0.2055, 0.2074, 0.2159, 0.2269, 0.2394,
0.2511, 0.2685, 0.2901, 0.3185, 0.3519, 0.3933, 0.4356, 0.4828, 0.5285,
0.5675, 0.5949, 0.6019, 0.5952, 0.5786, 0.5581, 0.5381, 0.5074, 0.4600,
0.4081, 0.3663, 0.3398, 0.3294, 0.3158, 0.3062, 0.2960, 0.2843, 0.2737,
0.2639, 0.2553, 0.2474, 0.2474};
static
float lu1_50[76] ={-.1746, -.1561, -.0735, 0.0810, 0.2161, 0.2430, 0.2246, 0.1953, 0.1629,
0.1355, 0.1127, 0.1027, 0.1005, 0.1007, 0.0977, 0.0944, 0.0899, 0.0874,
0.0906, 0.1065, 0.1776, 0.2971, 0.4673, 0.5043, 0.5361, 0.5538, 0.5589,
0.5556, 0.5432, 0.5300, 0.5190, 0.5068, 0.5010, 0.5010, 0.5045, 0.5090,
0.5140, 0.5187, 0.5235, 0.5427, 0.5573, 0.5663, 0.5696, 0.5665, 0.5601,
0.5470, 0.5313, 0.5112, 0.4848, 0.4526, 0.4145, 0.3718, 0.3248, 0.2765,
0.2280, 0.1866, 0.1539, 0.1375, 0.1340, 0.1392, 0.1484, 0.1622, 0.1831,
0.2186, 0.2632, 0.3045, 0.3382, 0.3621, 0.3798, 0.3936, 0.4060, 0.4195,
0.4319, 0.4433, 0.4538, 0.4633};
static
float lu2_50[76] ={0.7388, 0.7508, 0.7067, 0.5890, 0.4687, 0.4464, 0.4702, 0.5045, 0.5391,
0.5655, 0.5855, 0.5914, 0.5893, 0.5895, 0.5950, 0.6017, 0.6098, 0.6167,
0.6165, 0.5995, 0.5321, 0.3983, 0.1961, 0.1582, 0.1309, 0.1239, 0.1326,
0.1498, 0.1760, 0.2012, 0.2216, 0.2421, 0.2533, 0.2567, 0.2554, 0.2529,
0.2498, 0.2471, 0.2442, 0.2320, 0.2252, 0.2253, 0.2314, 0.2435, 0.2577,
0.2781, 0.2997, 0.3251, 0.3568, 0.3941, 0.4367, 0.4829, 0.5320, 0.5802,
0.6257, 0.6592, 0.6794, 0.6789, 0.6636, 0.6402, 0.6161, 0.5888, 0.5546,
0.5046, 0.4475, 0.3992, 0.3636, 0.3415, 0.3269, 0.3161, 0.3059, 0.2930,
0.2808, 0.2697, 0.2595, 0.2504};
SurfPtrS = Surface[0];
SurfPtrL = Surface[1];
if (Flags.first_pass == ON) {
/* create hash table for square root */
for (i = 0; i < (NSQR_TAB+1); ++i) {
sqtab[i] = 1.0 - (sqrt( (float)(i) / (float)NSQR_TAB ));
}
/* limb darkening coefficients */
for (i = 0; i < 2; ++i) {
Tn = Binary[i].Temperature;
if (Tn > SWITCH_TEMP) {
if (Binary[i].log_g <= 3.5) { lu1 = lu1_35; lu2 = lu2_35; }
else if (Binary[i].log_g == 4.0) { lu1 = lu1_40; lu2 = lu2_40; }
else if (Binary[i].log_g == 4.5) { lu1 = lu1_45; lu2 = lu2_45; }
else /* if (Binary[i].log_g >= 5.0) */ { lu1 = lu1_50; lu2 = lu2_50; }
k = 75; /* number of data in table */
lT = TK; /* pointer to temperature table */
} else {
if (Binary[i].log_g <= 3.5) { lu1 = hu1_35; lu2 = hu2_35; }
else if (Binary[i].log_g == 4.0) { lu1 = hu1_40; lu2 = hu2_40; }
else if (Binary[i].log_g == 4.5) { lu1 = hu1_45; lu2 = hu2_45; }
else /* if (Binary[i].log_g >= 5.0) */ { lu1 = hu1_50; lu2 = hu2_50; }
k = 34; /* number of data in table */
lT = TH; /* pointer to temperature table */
}
if (Tn <= lT[0]) {
u1[i] = lu1[0]; u2[i] = lu2[0];
}
else if (Tn >= lT[k]) {
u1[i] = lu1[k]; u2[i] = lu2[k];
}
else {
for (j = 0; j < k; ++j) {
if ( Tn >= lT[j] && Tn <= (lT[j+1]-DBL_EPSILON) ) {
u1[i] = lu1[j]; u2[i] = lu2[j];
break;
}
}
}
}
/* normalization for limb darkening coefficients */
NormLimb_[Primary] = 1.0/(PI*(1.0- u1[Primary]/3.0 - u2[Primary]/5.0));
NormLimb_[Secondary] = 1.0/(PI*(1.0- u1[Secondary]/3.0 - u2[Secondary]/5.0));
} /* end first_pass */
/* arrange properly for simple geometric test */
if (Binary[Primary].Radius >= Binary[Secondary].Radius) {
LargeStar = Primary;
SmallStar = Secondary;
} else {
LargeStar = Secondary;
SmallStar = Primary;
}
NumLarge = Binary[LargeStar].NumElem;
NumSmall = Binary[SmallStar].NumElem;
/* SmallStar is mirrored in y-z , ie. x -> (1-x) */
/* invisible side reached after Nmax non-sightings of other star */
Nmax = (4*StepsPerPi);
SurfPtrS = Surface[0];
SurfPtrL = Surface[1];
for (i = 0; i < MAXELEMENTS; ++i) {
Temp_[0][i] = SurfPtrS->temp; ++SurfPtrS;
Temp_[1][i] = SurfPtrL->temp; ++SurfPtrL;
}
for (k = 0; k < Flags.reflect; ++k) { /* multiple reflections */
#ifdef _WITH_GTK
if (Flags.interactive == ON) {
sprintf(message, _("Reflection iteration %1ld"), k);
my_appbar_push(message);
}
#endif
/* initialize the array for irradiated fluxes */
Init_ptr = &SumDelta_[0][0];
for (i = 0; i < 2*MAXELEMENTS; ++i) {
*Init_ptr = 0.0;
++Init_ptr;
}
LoopLow = (int)(NumLarge/3);
LoopHigh = (int)NumLarge;
#if defined(_WITH_MPI_COARSE)
/* all processes get here automatically
*/
if (Flags.elliptic != ON && numprocs > 1)
{
for (i = 0; i < MAXELEMENTS; ++i)
RecvBuf[i] = 0.0;
gsize = (numprocs < 100) ? numprocs : 100;
stride = (LoopHigh-LoopLow)/gsize;
displs = malloc(numprocs * sizeof(int));
rcount = malloc(numprocs * sizeof(int));
sdispls = malloc(numprocs * sizeof(int));
srcount = malloc(numprocs * sizeof(int));
if (displs == NULL || rcount == NULL ||
sdispls == NULL || srcount == NULL)
nf_error(_(errmsg[0]));
for (i = 0; i < numprocs; ++i)
{
if (i < gsize)
{
displs[i] = LoopLow + i * stride;
if ((displs[i] + rcount[i]) > NumLarge)
rcount[i] = NumLarge - displs[i];
rcount[i] = stride;
if ((i == (gsize-1)) && (displs[i] + rcount[i] < NumLarge))
rcount[i] = NumLarge - displs[i];
}
else
{
displs[i] = 0;
rcount[i] = 0;
}
}
LoopLow = displs[myrank];
LoopHigh = LoopLow + rcount[myrank];
for (i = 0; i < numprocs; ++i)
{
sdispls[i] = displs[myrank];
srcount[i] = rcount[myrank];
}
}
if (Flags.debug[VERBOSE] == ON)
fprintf(stderr, "MPI: Process %d of %d on %s: Reflect %d to %d\n",
myrank, numprocs, processor_name, LoopLow, LoopHigh-1);
#endif
SurfPtrL = &Surface[LargeStar][LoopLow];
/* -------- outer loop --------- */
i = LoopLow;
while (i < LoopHigh && myrank < gsize) {
/* simple geometric test, do we look towards other star ? */
if (SurfPtrL->lambda >= FLT_EPSILON) { /* yes, we do */
j = NumSmall - 1; Nnow = 0;
SurfPtrS = &Surface[SmallStar][j];
/* ---- inner loop --------- */
while ( ( j >= 0 ) && ( Nnow <= Nmax ) ) {
dx = SurfPtrL->lambda
- (1.0 - SurfPtrS->lambda);
dy = SurfPtrL->mu
- SurfPtrS->mu;
dz = SurfPtrL->nu
- SurfPtrS->nu;
/* we only need the correct sign, thus we do not divide */
/* by the norm of (ri-rj) */
/* thus cosgamma is really cosgamma*distance here */
cosgamma_ji = ( (-SurfPtrS->l) * dx
+ SurfPtrS->m * dy
+ SurfPtrS->n * dz );
if (cosgamma_ji >= DBL_EPSILON) {
cosgamma_ij = (SurfPtrL->l * (-dx)
+ SurfPtrL->m * (-dy)
+ SurfPtrL->n * (-dz) );
if (cosgamma_ij >= DBL_EPSILON) {
/* --------- visible, proceed ---------------------------- */
SurfPtrS->AreaType = FRONTSIDE;
SurfPtrL->AreaType = FRONTSIDE;
Nnow = 0; /* reset Nnow */
/* we only need the square of the distance */
distance_ij = (dx*dx) + (dy*dy) + (dz*dz);
/* ummm .. need distance for cosgamma because of */
/* the limb darkening ... */
distance_ij = sqrt(distance_ij);
cosgamma_ji = cosgamma_ji/distance_ij;
cosgamma_ij = cosgamma_ij/distance_ij;
/* in case of elliptic orbit */
if (Flags.elliptic == ON)
distance_ij = distance_ij * Orbit.Dist;
Delta_F = (cosgamma_ji * cosgamma_ij)
/ SQR(distance_ij);
/* irradiation onto LargeStar */
Tsmall = SurfPtrS->temp;
Tsmall = (Tsmall * Tsmall); /* T^2 */
Tsmall = (Tsmall * Tsmall); /* T^2 * T^2 = T^4 */
Delta_Ii = Delta_F
* SurfPtrS->area
* (1.0 - u1[SmallStar]*(1.0-cosgamma_ji)
- u2[SmallStar]*
sqtab[(int)(NSQR_TAB * cosgamma_ji)])
* (Tsmall);
/* irradiation on SmallStar */
Tlarge = SurfPtrL->temp;
Tlarge = (Tlarge * Tlarge); /* T^2 */
Tlarge = (Tlarge * Tlarge); /* T^2 * T^2 = T^4 */
Delta_Ij = Delta_F
* SurfPtrL->area
* (1.0 - u1[LargeStar]*(1.0-cosgamma_ij)
- u2[LargeStar]*
sqtab[(int)(NSQR_TAB * cosgamma_ij)])
* (Tlarge);
SumDelta_[LargeStar][i] = SumDelta_[LargeStar][i] + Delta_Ii;
SumDelta_[SmallStar][j] = SumDelta_[SmallStar][j] + Delta_Ij;
} else {
++Nnow;
}
} else {
++Nnow;
}
--j;
--SurfPtrS;
} /* end inner loop over surface elements */
} /* end if condition */
++SurfPtrL;
++i;
} /* end outer loop over surface elements */
/* ---- T(new)**4 = T(old)**4 + A * F(irr) --------------- */
/* for (i = 0; i < Binary[LargeStar].NumElem; ++i) */
for (i = LoopLow; i < LoopHigh; ++i)
{
Ttmp = Temp_[LargeStar][i] * Temp_[LargeStar][i];
/* Surface[LargeStar][i].temp = */
SumDelta_[LargeStar][i] =
sqrt(sqrt(
/* ( Temp_[Primary][i] * Temp_[Primary][i]
* Temp_[Primary][i] * Temp_[Primary][i] ) */
Ttmp * Ttmp
+ NormLimb_[SmallStar] * Binary[LargeStar].Albedo
* SumDelta_[LargeStar][i]
));
}
#if defined(_WITH_MPI_COARSE)
if (Flags.elliptic != ON && numprocs > 1)
{
MPI_Alltoallv (&SumDelta_[LargeStar][0], srcount, sdispls, MPI_DOUBLE,
&SumDelta_[LargeStar][0], rcount, displs, MPI_DOUBLE,
MPI_COMM_WORLD);
MPI_Allreduce(SumDelta_[SmallStar], RecvBuf, NumSmall, MPI_DOUBLE,
MPI_SUM, MPI_COMM_WORLD);
for (i = 0; i < NumSmall; ++i)
SumDelta_[SmallStar][i] = RecvBuf[i];
}
#endif
for (i = (NumLarge/3); i < NumLarge; ++i)
Surface[LargeStar][i].temp = SumDelta_[LargeStar][i];
for (i = 0; i < Binary[SmallStar].NumElem; ++i)
{
/* printf("N %4ld T %8.2f ", i, Surface[Secondary][i].temp); */
Ttmp = Temp_[SmallStar][i] * Temp_[SmallStar][i];
Surface[SmallStar][i].temp =
sqrt(sqrt(
/* ( Temp_[Secondary][i] * Temp_[Secondary][i]
* Temp_[Secondary][i] * Temp_[Secondary][i] ) */
Ttmp * Ttmp
+ NormLimb_[LargeStar] * Binary[SmallStar].Albedo
* SumDelta_[SmallStar][i]
));
}
} /* end multiple reflections */
}
syntax highlighted by Code2HTML, v. 0.9.1