1 #ifndef ERF_ORB_COS_ZENITH_H
2 #define ERF_ORB_COS_ZENITH_H
29 static constexpr
real dayspy =
real(365.0);
30 static constexpr
real ve =
real( 80.5);
51 lambm = lambm0 + (calday - ve)*2.0*
PI/dayspy;
58 lamb = lambm + eccen*( 2.0*sinl + eccen*( 1.25*sin(2.0*lmm)
59 + eccen*((13.0/12.0)*std::sin(3.0*lmm) - 0.25*sinl)));
69 invrho = (1.0 + eccen*std::cos(lamb - mvelpp)) / (1.0 - eccen*eccen);
73 delta = std::asin(std::sin(obliqr)*std::sin(lamb));
109 std::vector<real> obamp = {-2462.2214466, -857.3232075, -629.3231835,
110 -414.2804924, -311.7632587, 308.9408604,
111 -162.5533601, -116.1077911, 101.1189923,
112 -67.6856209, 24.9079067, 22.5811241,
113 -21.1648355, -15.6549876, 15.3936813,
114 14.6660938, -11.7273029, 10.2742696,
115 6.4914588, 5.8539148, -5.4872205,
116 -5.4290191, 5.1609570, 5.0786314,
117 -4.0735782, 3.7227167, 3.3971932,
118 -2.8347004, -2.6550721, -2.5717867,
119 -2.4712188, 2.4625410, 2.2464112,
120 -2.0755511, -1.9713669, -1.8813061,
121 -1.8468785, 1.8186742, 1.7601888,
122 -1.5428851, 1.4738838, -1.4593669,
123 1.4192259, -1.1818980, 1.1756474,
124 -1.1316126, 1.0896928};
126 std::vector<real> obrate = {31.609974, 32.620504, 24.172203,
127 31.983787, 44.828336, 30.973257,
128 43.668246, 32.246691, 30.599444,
129 42.681324, 43.836462, 47.439436,
130 63.219948, 64.230478, 1.010530,
131 7.437771, 55.782177, 0.373813,
132 13.218362, 62.583231, 63.593761,
133 76.438310, 45.815258, 8.448301,
134 56.792707, 49.747842, 12.058272,
135 75.278220, 65.241008, 64.604291,
136 1.647247, 7.811584, 12.207832,
137 63.856665, 56.155990, 77.448840,
138 6.801054, 62.209418, 20.656133,
139 48.344406, 55.145460, 69.000539,
140 11.071350, 74.291298, 11.047742,
141 0.636717, 12.844549};
144 std::vector<real> obphas = {251.9025, 280.8325, 128.3057,
145 292.7252, 15.3747, 263.7951,
146 308.4258, 240.0099, 222.9725,
147 268.7809, 316.7998, 319.6024,
148 143.8050, 172.7351, 28.9300,
149 123.5968, 20.2082, 40.8226,
150 123.4722, 155.6977, 184.6277,
151 267.2772, 55.0196, 152.5268,
152 49.1382, 204.6609, 56.5233,
153 200.3284, 201.6651, 213.5577,
154 17.0374, 164.4194, 94.5422,
155 131.9124, 61.0309, 296.2073,
156 135.4894, 114.8750, 247.0691,
157 256.6114, 32.1008, 143.6804,
158 16.8784, 160.6835, 27.5932,
166 std::vector<real> ecamp = { 0.01860798, 0.01627522, -0.01300660,
167 0.00988829, -0.00336700, 0.00333077,
168 -0.00235400, 0.00140015, 0.00100700,
169 0.00085700, 0.00064990, 0.00059900,
170 0.00037800, -0.00033700, 0.00027600,
171 0.00018200, -0.00017400, -0.00012400,
175 std::vector<real> ecrate = { 4.2072050, 7.3460910, 17.8572630,
176 17.2205460, 16.8467330, 5.1990790,
177 18.2310760, 26.2167580, 6.3591690,
178 16.2100160, 3.0651810, 16.5838290,
179 18.4939800, 6.1909530, 18.8677930,
180 17.4255670, 6.1860010, 18.4174410,
184 std::vector<real> ecphas = { 28.620089, 193.788772, 308.307024,
185 320.199637, 279.376984, 87.195000,
186 349.129677, 128.443387, 154.143880,
187 291.269597, 114.860583, 332.092251,
188 296.414411, 145.769910, 337.237063,
189 152.092288, 126.839891, 210.667199,
196 std::vector<real> mvamp = { 7391.0225890, 2555.1526947, 2022.7629188,
197 -1973.6517951, 1240.2321818, 953.8679112,
198 -931.7537108, 872.3795383, 606.3544732,
199 -496.0274038, 456.9608039, 346.9462320,
200 -305.8412902, 249.6173246, -199.1027200,
201 191.0560889, -175.2936572, 165.9068833,
202 161.1285917, 139.7878093, -133.5228399,
203 117.0673811, 104.6907281, 95.3227476,
204 86.7824524, 86.0857729, 70.5893698,
205 -69.9719343, -62.5817473, 61.5450059,
206 -57.9364011, 57.1899832, -57.0236109,
207 -54.2119253, 53.2834147, 52.1223575,
208 -49.0059908, -48.3118757, -45.4191685,
209 -42.2357920, -34.7971099, 34.4623613,
210 -33.8356643, 33.6689362, -31.2521586,
211 -30.8798701, 28.4640769, -27.1960802,
212 27.0860736, -26.3437456, 24.7253740,
213 24.6732126, 24.4272733, 24.0127327,
214 21.7150294, -21.5375347, 18.1148363,
215 -16.9603104, -16.1765215, 15.5567653,
216 15.4846529, 15.2150632, 14.5047426,
217 -14.3873316, 13.1351419, 12.8776311,
218 11.9867234, 11.9385578, 11.7030822,
219 11.6018181, -11.2617293, -10.4664199,
220 10.4333970, -10.2377466, 10.1934446,
221 -10.1280191, 10.0289441, -10.0034259};
224 std::vector<real> mvrate = {31.609974, 32.620504, 24.172203,
225 0.636717, 31.983787, 3.138886,
226 30.973257, 44.828336, 0.991874,
227 0.373813, 43.668246, 32.246691,
228 30.599444, 2.147012, 10.511172,
229 42.681324, 13.650058, 0.986922,
230 9.874455, 13.013341, 0.262904,
231 0.004952, 1.142024, 63.219948,
232 0.205021, 2.151964, 64.230478,
233 43.836462, 47.439436, 1.384343,
234 7.437771, 18.829299, 9.500642,
235 0.431696, 1.160090, 55.782177,
236 12.639528, 1.155138, 0.168216,
237 1.647247, 10.884985, 5.610937,
238 12.658184, 1.010530, 1.983748,
239 14.023871, 0.560178, 1.273434,
240 12.021467, 62.583231, 63.593761,
241 76.438310, 4.280910, 13.218362,
242 17.818769, 8.359495, 56.792707,
243 8.448301, 1.978796, 8.863925,
244 0.186365, 8.996212, 6.771027,
245 45.815258, 12.002811, 75.278220,
246 65.241008, 18.870667, 22.009553,
247 64.604291, 11.498094, 0.578834,
248 9.237738, 49.747842, 2.147012,
249 1.196895, 2.133898, 0.173168};
252 std::vector<real> mvphas = {251.9025, 280.8325, 128.3057,
253 348.1074, 292.7252, 165.1686,
254 263.7951, 15.3747, 58.5749,
255 40.8226, 308.4258, 240.0099,
256 222.9725, 106.5937, 114.5182,
257 268.7809, 279.6869, 39.6448,
258 126.4108, 291.5795, 307.2848,
259 18.9300, 273.7596, 143.8050,
260 191.8927, 125.5237, 172.7351,
261 316.7998, 319.6024, 69.7526,
262 123.5968, 217.6432, 85.5882,
263 156.2147, 66.9489, 20.2082,
264 250.7568, 48.0188, 8.3739,
265 17.0374, 155.3409, 94.1709,
266 221.1120, 28.9300, 117.1498,
267 320.5095, 262.3602, 336.2148,
268 233.0046, 155.6977, 184.6277,
269 267.2772, 78.9281, 123.4722,
270 188.7132, 180.1364, 49.1382,
271 152.5268, 98.2198, 97.4808,
272 221.5376, 168.2438, 161.1199,
273 55.0196, 262.6495, 200.3284,
274 201.6651, 294.6547, 99.8233,
275 213.5577, 154.1631, 232.7153,
276 138.3034, 204.6609, 106.5938,
277 250.4676, 332.3345, 27.3039};
298 eccen2 = eccen*eccen;
299 eccen3 = eccen2*eccen;
328 yb4_1950AD =
real(1950.0) -
real(iyear_AD);
329 years = - yb4_1950AD;
343 for (
int i(0); i<poblen; ++i) {
344 obsum = obsum + obamp[i]*psecdeg*std::cos( (obrate[i]*psecdeg*years + obphas[i]) *
degrad );
346 obliq =
real(23.320556) + obsum;
357 for (
int i(0); i<pecclen; ++i) {
358 cossum = cossum + ecamp[i]*std::cos( (ecrate[i]*psecdeg*years+ecphas[i]) *
degrad );
362 for (
int i(0); i<pecclen; ++i) {
363 sinsum = sinsum + ecamp[i]*std::sin( (ecrate[i]*psecdeg*years+ecphas[i]) *
degrad );
368 eccen2 = cossum*cossum + sinsum*sinsum;
369 eccen = std::sqrt(eccen2);
370 eccen3 = eccen2*eccen;
373 if (std::fabs(cossum) <= 1.0e-8) {
376 }
else if (sinsum <= 0.0) {
378 }
else if (sinsum > 0.0) {
381 }
else if (cossum <= 0.0) {
382 fvelp = std::atan(sinsum/cossum) +
PI;
385 fvelp = std::atan(sinsum/cossum) + 2.0*
PI;
387 fvelp = std::atan(sinsum/cossum);
402 for (
int i(0); i<pmvelen; ++i) {
403 mvsum = mvsum + mvamp[i]*psecdeg*std::sin( (mvrate[i]*psecdeg*years + mvphas[i]) *
degrad);
405 mvelp = fvelp/
degrad +
real(50.439273)*psecdeg*years +
real(3.392506) + mvsum;
410 }
while (mvelp < 0.0);
414 }
while (mvelp >= 360.0);
434 mvelpp = (mvelp + 180.0)*
degrad;
438 beta = std::sqrt(1.0 - eccen2);
447 lambm0 = 2.0*( (.5*eccen + .125*eccen3)*(1.0 + beta)*std::sin(mvelpp)
448 - .250*eccen2*(.5 + beta)*std::sin(2.0*mvelpp)
449 + .125*eccen3*(1./3. + beta)*std::sin(3.0*mvelpp) );
465 del = lat -
real(1.0e-05);
466 }
else if (lat == -
PIoTwo) {
467 del = lat +
real(1.0e-05);
475 phi = declin -
real(1.0e-05);
476 }
else if (declin == -
PIoTwo) {
477 phi = declin +
real(1.0e-05);
485 real cos_h = - std::tan(del) * std::tan(phi);
488 }
else if (cos_h >= 1.0) {
491 h = std::acos(cos_h);
496 real t1 = (jday - int(jday)) * 2.0*
PI + lon -
PI;
499 }
else if (t1 < -
PI) {
508 real aa = std::sin(lat) * std::sin(declin);
509 real bb = std::cos(lat) * std::cos(declin);
514 real tt1,tt2,tt3,tt4;
515 if ( (t2 >=
PI) && (t1 <=
PI) && (
PI - h <= dt) ) {
517 tt1 = std::min(std::max(t1, -h), h);
518 tt4 = std::min(std::max(t2, 2.0*
PI - h), 2.0*
PI + h);
520 }
else if ( (t2 >= -
PI) && (t1 <= -
PI) && (
PI - h <= dt) ) {
522 tt1 = std::min(std::max(t1, -2.0*
PI - h), -2.0*
PI + h);
523 tt4 = std::min(std::max(t2, -h), h);
527 tt2 = std::min(std::max(t2 - 2.0*
PI, -h), h);
528 }
else if (t2 < -
PI) {
529 tt2 = std::min(std::max(t2 + 2.0*
PI, -h), h);
531 tt2 = std::min(std::max(t2 , -h), h);
535 tt1 = std::min(std::max(t1 - 2.0*
PI, -h), h);
536 }
else if (t1 < -
PI) {
537 tt1 = std::min(std::max(t1 + 2.0*
PI, -h), h);
539 tt1 = std::min(std::max(t1 , -h), h);
547 if ( (tt2 > tt1) || (tt4 > tt3) ) {
548 return (aa * (tt2 - tt1) + bb * (sin(tt2) - sin(tt1))) / dt +
549 (aa * (tt4 - tt3) + bb * (sin(tt4) - sin(tt3))) / dt;
564 real uniform_angle = -1.,
565 real constant_zenith_angle_deg = -1.)
568 if ( constant_zenith_angle_deg >= 0. ) {
569 return std::cos( constant_zenith_angle_deg *
PI/180. );
573 if ( uniform_angle >= 0.) {
574 return std::cos(uniform_angle);
578 bool use_dt_avg =
false;
586 return std::sin(lat)*std::sin(declin) - std::cos(lat)*std::cos(declin) *
587 std::cos((jday-floor(jday))*
real(2.0)*
PI + lon);
static constexpr int ORB_UNDEF_INT
Definition: ERF_Constants.H:96
constexpr amrex::Real PI
Definition: ERF_Constants.H:6
constexpr amrex::Real PIoTwo
Definition: ERF_Constants.H:7
double real
Definition: ERF_OrbCosZenith.H:9
AMREX_GPU_HOST AMREX_FORCE_INLINE real orbital_cos_zenith(real &jday, real &lat, real &lon, real &declin, real dt_avg=-1., real uniform_angle=-1., real constant_zenith_angle_deg=-1.)
Definition: ERF_OrbCosZenith.H:559
AMREX_GPU_HOST AMREX_FORCE_INLINE real orbital_avg_cos_zenith(real &jday, real &lat, real &lon, real &declin, real &dt_avg)
Definition: ERF_OrbCosZenith.H:456
AMREX_GPU_HOST AMREX_FORCE_INLINE void orbital_decl(real &calday, real &eccen, real &mvelpp, real &lambm0, real &obliqr, real &delta, real &eccf)
Definition: ERF_OrbCosZenith.H:15
AMREX_GPU_HOST AMREX_FORCE_INLINE void orbital_params(int &iyear_AD, real &eccen, real &obliq, real &mvelp, real &obliqr, real &lambm0, real &mvelpp)
Definition: ERF_OrbCosZenith.H:81
real(c_double), parameter degrad
Definition: ERF_module_model_constants.F90:75