1/** 2 * Rise and Set of Sun and Moon of Earth planet in current galaxy 3 * (address: Local Bubble in Local Interstellar Cloud of OrionâCygnus Arm, for more coordinats ask Frank the Pug or agent K) 4 */ 5 6no_rs = { 7 hours: '--', 8 minutes: '--' 9 }; 10 11function GetRiseSet(glat,glong) { 12 var OutString = ""; 13 var calend; 14 var quady = new Array; 15 var sunp = new Array; 16 var moonp = new Array; 17 var y, m, day, tz, mj, lst1, i; 18 var rads = 0.0174532925, sinmoonalt; 19 // 20 // parse the form to make sure the numbers are numbers and not strings! 21 // 22 var d = new Date(); 23 y = d.getFullYear(); 24 m = d.getMonth()+1; 25 day = d.getDate(); 26 tz = -d.getTimezoneOffset()/60; 27 28 // 29 // main loop. All the work is done in the functions with the long names 30 // find_sun_and_twi_events_for_date() and find_moonrise_set() 31 // 32 mj = mjd(day, m, y, 0.0); 33 34 return { 35 sun: find_sun_and_twi_events_for_date(mj, tz, glong, glat), 36 moon: find_moonrise_set(mj, tz, glong, glat) 37 }; 38 } // end of main program 39 40 41function hrsmin(hours) { 42// 43// takes decimal hours and returns a string in hhmm format 44// 45 var hrs, h, m, dum; 46 hrs = Math.floor(hours * 60 + 0.5)/ 60.0; 47 h = Math.floor(hrs); 48 m = Math.floor(60 * (hrs - h) + 0.5); 49 dum = String(h) + ':' + String(m); 50 // 51 // the jiggery pokery below is to make sure that two minutes past midnight 52 // comes out as 0002 not 2. Javascript does not appear to have 'format codes' 53 // like C 54 // 55 56 if(m < 10) 57 m = '0'+String(m); 58 if(h < 10) 59 h = '0'+String(h); 60 61 return { 62 hours: h, 63 minutes: m 64 }; 65 } 66 67 68function ipart(x) { 69// 70// returns the integer part - like int() in basic 71// 72 var a; 73 if (x> 0) { 74 a = Math.floor(x); 75 } 76 else { 77 a = Math.ceil(x); 78 } 79 return a; 80 } 81 82 83function frac(x) { 84// 85// returns the fractional part of x as used in minimoon and minisun 86// 87 var a; 88 a = x - Math.floor(x); 89 if (a < 0) a += 1; 90 return a; 91 } 92 93// 94// round rounds the number num to dp decimal places 95// the second line is some C like jiggery pokery I 96// found in an OReilly book which means if dp is null 97// you get 2 decimal places. 98// 99 function round(num, dp) { 100// dp = (!dp ? 2: dp); 101 return Math.round (num * Math.pow(10, dp)) / Math.pow(10, dp); 102 } 103 104 105function range(x) { 106// 107// returns an angle in degrees in the range 0 to 360 108// 109 var a, b; 110 b = x / 360; 111 a = 360 * (b - ipart(b)); 112 if (a < 0 ) { 113 a = a + 360 114 } 115 return a 116 } 117 118 119function mjd(day, month, year, hour) { 120// 121// Takes the day, month, year and hours in the day and returns the 122// modified julian day number defined as mjd = jd - 2400000.5 123// checked OK for Greg era dates - 26th Dec 02 124// 125 var a, b; 126 if (month <= 2) { 127 month = month + 12; 128 year = year - 1; 129 } 130 a = 10000.0 * year + 100.0 * month + day; 131 if (a <= 15821004.1) { 132 b = -2 * Math.floor((year + 4716)/4) - 1179; 133 } 134 else { 135 b = Math.floor(year/400) - Math.floor(year/100) + Math.floor(year/4); 136 } 137 a = 365.0 * year - 679004.0; 138 return (a + b + Math.floor(30.6001 * (month + 1)) + day + hour/24.0); 139 } 140 141function caldat(mjd) { 142// 143// Takes mjd and returns the civil calendar date in Gregorian calendar 144// as a string in format yyyymmdd.hhhh 145// looks OK for Greg era dates - not good for earlier - 26th Dec 02 146// 147 var calout; 148 var b, d, f, jd, jd0, c, e, day, month, year, hour; 149 jd = mjd + 2400000.5; 150 jd0 = Math.floor(jd + 0.5); 151 if (jd0 < 2299161.0) { 152 c = jd0 + 1524.0; 153 } 154 else { 155 b = Math.floor((jd0 - 1867216.25) / 36524.25); 156 c = jd0 + (b - Math.floor(b/4)) + 1525.0; 157 } 158 d = Math.floor((c - 122.1)/365.25); 159 e = 365.0 * d + Math.floor(d/4); 160 f = Math.floor(( c - e) / 30.6001); 161 day = Math.floor(c - e + 0.5) - Math.floor(30.6001 * f); 162 month = f - 1 - 12 * Math.floor(f/14); 163 year = d - 4715 - Math.floor((7 + month)/10); 164 hour = 24.0 * (jd + 0.5 - jd0); 165 hour = hrsmin(hour); 166 calout = round(year * 10000.0 + month * 100.0 + day + hour/10000, 4); 167 return calout + ""; //making sure calout is a string 168 } 169 170 171function quad(ym, yz, yp) { 172// 173// finds the parabola throuh the three points (-1,ym), (0,yz), (1, yp) 174// and returns the coordinates of the max/min (if any) xe, ye 175// the values of x where the parabola crosses zero (roots of the quadratic) 176// and the number of roots (0, 1 or 2) within the interval [-1, 1] 177// 178// well, this routine is producing sensible answers 179// 180// results passed as array [nz, z1, z2, xe, ye] 181// 182 var nz, a, b, c, dis, dx, xe, ye, z1, z2, nz; 183 var quadout = new Array; 184 185 nz = 0; 186 a = 0.5 * (ym + yp) - yz; 187 b = 0.5 * (yp - ym); 188 c = yz; 189 xe = -b / (2 * a); 190 ye = (a * xe + b) * xe + c; 191 dis = b * b - 4.0 * a * c; 192 if (dis > 0) { 193 dx = 0.5 * Math.sqrt(dis) / Math.abs(a); 194 z1 = xe - dx; 195 z2 = xe + dx; 196 if (Math.abs(z1) <= 1.0) nz += 1; 197 if (Math.abs(z2) <= 1.0) nz += 1; 198 if (z1 < -1.0) z1 = z2; 199 } 200 quadout[0] = nz; 201 quadout[1] = z1; 202 quadout[2] = z2; 203 quadout[3] = xe; 204 quadout[4] = ye; 205 return quadout; 206 } 207 208 209function lmst(mjd, glong) { 210// 211// Takes the mjd and the longitude (west negative) and then returns 212// the local sidereal time in hours. Im using Meeus formula 11.4 213// instead of messing about with UTo and so on 214// 215 var lst, t, d; 216 d = mjd - 51544.5 217 t = d / 36525.0; 218 lst = range(280.46061837 + 360.98564736629 * d + 0.000387933 *t*t - t*t*t / 38710000); 219 return (lst/15.0 + glong/15); 220 } 221 222 223function minisun(t) { 224// 225// returns the ra and dec of the Sun in an array called suneq[] 226// in decimal hours, degs referred to the equinox of date and using 227// obliquity of the ecliptic at J2000.0 (small error for +- 100 yrs) 228// takes t centuries since J2000.0. Claimed good to 1 arcmin 229// 230 var p2 = 6.283185307, coseps = 0.91748, sineps = 0.39778; 231 var L, M, DL, SL, X, Y, Z, RHO, ra, dec; 232 var suneq = new Array; 233 234 M = p2 * frac(0.993133 + 99.997361 * t); 235 DL = 6893.0 * Math.sin(M) + 72.0 * Math.sin(2 * M); 236 L = p2 * frac(0.7859453 + M / p2 + (6191.2 * t + DL)/1296000); 237 SL = Math.sin(L); 238 X = Math.cos(L); 239 Y = coseps * SL; 240 Z = sineps * SL; 241 RHO = Math.sqrt(1 - Z * Z); 242 dec = (360.0 / p2) * Math.atan(Z / RHO); 243 ra = (48.0 / p2) * Math.atan(Y / (X + RHO)); 244 if (ra <0 ) ra += 24; 245 suneq[1] = dec; 246 suneq[2] = ra; 247 return suneq; 248 } 249 250 251function minimoon(t) { 252// 253// takes t and returns the geocentric ra and dec in an array mooneq 254// claimed good to 5' (angle) in ra and 1' in dec 255// tallies with another approximate method and with ICE for a couple of dates 256// 257 var p2 = 6.283185307, arc = 206264.8062, coseps = 0.91748, sineps = 0.39778; 258 var L0, L, LS, F, D, H, S, N, DL, CB, L_moon, B_moon, V, W, X, Y, Z, RHO; 259 var mooneq = new Array; 260 261 L0 = frac(0.606433 + 1336.855225 * t); // mean longitude of moon 262 L = p2 * frac(0.374897 + 1325.552410 * t) //mean anomaly of Moon 263 LS = p2 * frac(0.993133 + 99.997361 * t); //mean anomaly of Sun 264 D = p2 * frac(0.827361 + 1236.853086 * t); //difference in longitude of moon and sun 265 F = p2 * frac(0.259086 + 1342.227825 * t); //mean argument of latitude 266 267 // corrections to mean longitude in arcsec
268 DL = 22640 * Math.sin(L) 269 DL += -4586 * Math.sin(L - 2*D); 270 DL += +2370 * Math.sin(2*D); 271 DL += +769 * Math.sin(2*L); 272 DL += -668 * Math.sin(LS); 273 DL += -412 * Math.sin(2*F); 274 DL += -212 * Math.sin(2*L - 2*D); 275 DL += -206 * Math.sin(L + LS - 2*D); 276 DL += +192 * Math.sin(L + 2*D); 277 DL += -165 * Math.sin(LS - 2*D); 278 DL += -125 * Math.sin(D); 279 DL += -110 * Math.sin(L + LS); 280 DL += +148 * Math.sin(L - LS); 281 DL += -55 * Math.sin(2*F - 2*D); 282 283 // simplified form of the latitude terms 284 S = F + (DL + 412 * Math.sin(2*F) + 541* Math.sin(LS)) / arc; 285 H = F - 2*D; 286 N = -526 * Math.sin(H); 287 N += +44 * Math.sin(L + H); 288 N += -31 * Math.sin(-L + H); 289 N += -23 * Math.sin(LS + H); 290 N += +11 * Math.sin(-LS + H); 291 N += -25 * Math.sin(-2*L + F); 292 N += +21 * Math.sin(-L + F); 293 294 // ecliptic long and lat of Moon in rads 295 L_moon = p2 * frac(L0 + DL / 1296000); 296 B_moon = (18520.0 * Math.sin(S) + N) /arc; 297 298 // equatorial coord conversion - note fixed obliquity 299 CB = Math.cos(B_moon); 300 X = CB * Math.cos(L_moon); 301 V = CB * Math.sin(L_moon); 302 W = Math.sin(B_moon); 303 Y = coseps * V - sineps * W; 304 Z = sineps * V + coseps * W 305 RHO = Math.sqrt(1.0 - Z*Z); 306 dec = (360.0 / p2) * Math.atan(Z / RHO); 307 ra = (48.0 / p2) * Math.atan(Y / (X + RHO)); 308 if (ra <0 ) ra += 24; 309 mooneq[1] = dec; 310 mooneq[2] = ra; 311 return mooneq; 312 } 313 314 315function sin_alt(iobj, mjd0, hour, glong, cglat, sglat) { 316// 317// this rather mickey mouse function takes a lot of 318// arguments and then returns the sine of the altitude of 319// the object labelled by iobj. iobj = 1 is moon, iobj = 2 is sun 320// 321 var mjd, t, ra, dec, tau, salt, rads = 0.0174532925; 322 var objpos = new Array; 323 mjd = mjd0 + hour/24.0; 324 t = (mjd - 51544.5) / 36525.0; 325 if (iobj == 1) { 326 objpos = minimoon(t); 327 } 328 else { 329 objpos = minisun(t); 330 } 331 ra = objpos[2]; 332 dec = objpos[1]; 333 // hour angle of object 334 tau = 15.0 * (lmst(mjd, glong) - ra); 335 // sin(alt) of object using the conversion formulas 336 salt = sglat * Math.sin(rads*dec) + cglat * Math.cos(rads*dec) * Math.cos(rads*tau); 337 return salt; 338 } 339 340 341function find_sun_and_twi_events_for_date(mjd, tz, glong, glat) { 342// 343// this is my attempt to encapsulate most of the program in a function 344// then this function can be generalised to find all the Sun events. 345// 346// 347 var sglong, sglat, date, ym, yz, above, utrise, utset, j; 348 var yp, nz, rise, sett, hour, z1, z2, iobj, rads = 0.0174532925; 349 var quadout = new Array; 350 var sinho = new Array; 351 var always_up = " ****"; 352 var always_down = " ...."; 353 var outstring = ""; 354// 355// Set up the array with the 4 values of sinho needed for the 4 356// kinds of sun event 357// 358 sinho[0] = Math.sin(rads * -0.833); //sunset upper limb simple refraction 359 sinho[1] = Math.sin(rads * -6.0); //civil twi 360 sinho[2] = Math.sin(rads * -12.0); //nautical twi 361 sinho[3] = Math.sin(rads * -18.0); //astro twi 362 sglat = Math.sin(rads * glat); 363 cglat = Math.cos(rads * glat); 364 date = mjd - tz/24; 365// 366// main loop takes each value of sinho in turn and finds the rise/set 367// events associated with that altitude of the Sun 368// 369 j = 0; 370 rise = false; 371 sett = false; 372 above = false; 373 hour = 1.0; 374 ym = sin_alt(2, date, hour - 1.0, glong, cglat, sglat) - sinho[j]; 375 if (ym > 0.0) above = true; 376 // 377 // the while loop finds the sin(alt) for sets of three consecutive 378 // hours, and then tests for a single zero crossing in the interval 379 // or for two zero crossings in an interval or for a grazing event 380 // The flags rise and sett are set accordingly 381 // 382 while(hour < 25 && (sett == false || rise == false)) { 383 yz = sin_alt(2, date, hour, glong, cglat, sglat) - sinho[j]; 384 yp = sin_alt(2, date, hour + 1.0, glong, cglat, sglat) - sinho[j]; 385 quadout = quad(ym, yz, yp); 386 nz = quadout[0]; 387 z1 = quadout[1]; 388 z2 = quadout[2]; 389 xe = quadout[3]; 390 ye = quadout[4]; 391 392 // case when one event is found in the interval 393 if (nz == 1) { 394 if (ym < 0.0) { 395 utrise = hour + z1; 396 rise = true; 397 } 398 else { 399 utset = hour + z1; 400 sett = true; 401 } 402 } // end of nz = 1 case 403
404 // case where two events are found in this interval 405 // (rare but whole reason we are not using simple iteration) 406 if (nz == 2) { 407 if (ye < 0.0) { 408 utrise = hour + z2; 409 utset = hour + z1; 410 } 411 else { 412 utrise = hour + z1; 413 utset = hour + z2; 414 } 415 } // end of nz = 2 case 416 417 // set up the next search interval 418 ym = yp; 419 hour += 2.0; 420 421 } // end of while loop 422 // 423 // now search has completed, we compile the string to pass back 424 // to the main loop. The string depends on several combinations 425 // of the above flag (always above or always below) and the rise 426 // and sett flags 427 // 428 out= { 429 rise: no_rs, 430 set: no_rs 431 }; 432 433 if (rise == true) out.rise = hrsmin(utrise); 434 if (sett == true) out.set = hrsmin(utset); 435 436 return out; 437 } 438 439function find_moonrise_set(mjd, tz, glong, glat) { 440// 441// Im using a separate function for moonrise/set to allow for different tabulations 442// of moonrise and sun events ie weekly for sun and daily for moon. The logic of 443// the function is identical to find_sun_and_twi_events_for_date() 444// 445 var sglong, sglat, date, ym, yz, above, utrise, utset, j; 446 var yp, nz, rise, sett, hour, z1, z2, iobj, rads = 0.0174532925; 447 var quadout = new Array; 448 var sinho; 449 var always_up = " ****"; 450 var always_down = " ...."; 451 var outstring = ""; 452 453 sinho = Math.sin(rads * 8/60); //moonrise taken as centre of moon at +8 arcmin 454 sglat = Math.sin(rads * glat); 455 cglat = Math.cos(rads * glat); 456 date = mjd - tz/24; 457 rise = false; 458 sett = false; 459 above = false; 460 hour = 1.0; 461 ym = sin_alt(1, date, hour - 1.0, glong, cglat, sglat) - sinho; 462 if (ym > 0.0) above = true; 463 while(hour < 25 && (sett == false || rise == false)) { 464 yz = sin_alt(1, date, hour, glong, cglat, sglat) - sinho; 465 yp = sin_alt(1, date, hour + 1.0, glong, cglat, sglat) - sinho; 466 quadout = quad(ym, yz, yp); 467 nz = quadout[0]; 468 z1 = quadout[1]; 469 z2 = quadout[2]; 470 xe = quadout[3]; 471 ye = quadout[4]; 472 473 // case when one event is found in the interval 474 if (nz == 1) { 475 if (ym < 0.0) { 476 utrise = hour + z1; 477 rise = true; 478 } 479 else { 480 utset = hour + z1; 481 sett = true; 482 } 483 } // end of nz = 1 case 484
485 // case where two events are found in this interval 486 // (rare but whole reason we are not using simple iteration) 487 if (nz == 2) { 488 if (ye < 0.0) { 489 utrise = hour + z2; 490 utset = hour + z1; 491 } 492 else { 493 utrise = hour + z1; 494 utset = hour + z2; 495 } 496 } 497 498 // set up the next search interval 499 ym = yp; 500 hour += 2.0; 501 502 } // end of while loop 503 504 out= { 505 rise: no_rs, 506 set: no_rs 507 }; 508 509 if (rise == true) out.rise = hrsmin(utrise); 510 if (sett == true) out.set = hrsmin(utset); 511 512 return out; 513}
Line numbers count LF bytes from the start of the resource, as the search results do. Vendor segments are library code the classifier recognised; they are stored but not indexed. Bytes are shown as Latin1 characters, one per byte.