/* 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" /**************************************************************************** @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 */ }