****************************************************************************** program testime ****************************************************************************** * tests some of the time and coordinate conversion routines * must be linked with timesub.for and update.for * implicit real*8 (a-h, o-z) real*8 l, month, hour, min, jd, lat, long, lst, lam character*72 msg, blank character*1 blanks(72) equivalence (blanks(1), blank) data blanks/72*' '/ type *,' Program Test Time Checks coordinate transforms' type *,' and will permit hand calculations of coordinates' type *,' ' data alp_src/5.5972/,del_src/-69.29933/ alp = alp_src del = del_src data detlat/36.42561/,detlong/ -137.3103/,zone/9./ c default coords are for superk location iopt = 1 10 continue type *, ' ' type *,' select option:' type *,' 0 = quit' type *,' 1 = equatorial to galactic coordinates' type *,' 2 = galactic to equatorial coordinates' type *,' 3 = conversion to julian date' type *,' 4 = conversion from date' type *,' 5 = local time to local siderial time' type *,' 6 = local siderial time to local time' type *,' 7 = local detector coords to equatorial coords' type *,' 8 = equatorial to local detector coords' type *,' 9 = equatorial to ecliptic coordinates' type *,' 10 = ecliptic to equatorial coordinates' type *,' 11 = local to solar coordinates' type *,' 12 = solar to local coordinates' type *,' 13 = get lunar ra and dec' type *,' 14 = get angle from moon' type *,' 15 = set detector and source coords' type *,' 16 = sun ra & dec and barycentric correction' type *,' 17 = cnvt evt coords to ra & dec, calc src ang dist' msg = ' option please? $' call iupdate(iopt,msg) goto(100, 200, 300, 400, 500, 600, 700, 800, 900, & 1000, 1100, 1200, 1300, 1400, 1500, 1600, 1700) iopt * otherwise quit stop 100 continue msg = ' ra in decimal hours $' call update(alp,msg) msg = 'declination in decimal deg $' call update(del,msg) call eqgal(alp,del,gl,gb) type *,'galactic longitude =',gl type *,'galactic latitude =',gb go to 10 200 continue msg = 'gal long in deg $' call update(gl,msg) msg = 'gal lat in deg $' call update(gb,msg) call galeq(gl,gb,alp,del) type *,' right ascension = ',alp,' hours' type *,' declination = ',del,' degrees' go to 10 300 continue msg = 'Greenwich day of month, 1-31$' call update(gday,msg) msg = 'G month, 1-12$' call update(gmonth,msg) msg = 'G year, ad$' call update(gyear,msg) call dayjul(gday,gmonth,gyear,rjd) type 350, rjd 350 format(' Julian day at 0 hrs GMT= ',f15.5) jd = rjd go to 10 400 continue msg = ' decimal Julian day = $' call update(jd,msg) call julday(jd,gyear,gmonth,gday,dgmt) type *,' G year = ',gyear type *,' G month = ',gmonth type *,' G day = ',gday type *,' GMT, hrs = ',dgmt go to 10 500 continue * here to convert local time to local siderial time type *,' enter detector time zone, allowing for dst' type *,' negative for West of Greenwich' msg = ' zone = $' call update(zone,msg) msg = ' det lat (deg N+) = $' call update(detlat,msg) msg = ' det long (deg W+) = $' call update(detlong,msg) msg = ' local year = $' call update(year,msg) msg = ' local month = $' call update(month,msg) msg = ' local day = $' call update(day,msg) msg = ' local hour = $' call update(hour,msg) msg = ' local minutes = $' call update(min,msg) msg = ' local seconds = $' call update(sec,msg) call ltlst(year, month, day, hour, min, sec, & detlat, detlong, zone, dlst) type *,' local sidereal time is = ',dlst,' hours' go to 10 600 continue * do conversion from local sidereal back to local time type *,' Enter detector time zone, allowing for dst,' type *,' negative for West of Greenwich.' msg = ' zone = $' call update(zone,msg) msg = ' local year = $' call update(year,msg) msg = ' local month = $' call update(month,msg) msg = ' local day = $' call update(day,msg) msg = ' local sidereal time, hrs = $' call update(dlst,msg) call lstlt(dlst, detlat, detlong, zone, & year, month, day, hour, min, sec) type *,' local time is = ',hour,' hours, ',min,' min, ',sec,' sec' go to 10 700 continue * convert local detector coords to equatorial coords msg = 'zenith angle re detector, deg$' call update(theta, msg) msg = 'azimuth re det, deg$' call update(phi,msg) msg = 'local sid time, decimal hrs$' call update(dlst,msg) msg = ' det lat = $' call update(detlat,msg) call loceq(theta, phi, dlst, detlat, alp, del) type *,' right ascension = ',alp,' hours' type *,' declination = ',del,' degrees' go to 10 800 continue * convert equatorial to local coordinates msg = 'ra, decimal hrs$' call update(alp, msg) msg = 'declination, deg$' call update(del, msg) msg = 'loc sid time, hrs$' call update(dlst, msg) msg = ' det lat = $' call update(detlat, msg) call eqloc(alp, del, dlst, detlat, theta, phi) type *,' zenith angle = ',theta,' degrees' type *,' azimuth, n of e =',phi,' degrees' go to 10 900 continue ! here for eq -> ecl msg = 'right ascension, hrs$' call update(alp, msg) msg = 'declination, degrees$' call update(del, msg) msg = 'year, AD$' call update(year, msg) call eqecl(alp, del, year, el, eb) type *,' ecliptic longitude = ',el,' degrees' type *,' ecliptic latitude = ',eb,' degrees' go to 10 1000 continue ! here for ecl -> eq msg = 'ecliptic latitude, degrees = $' call update(el, msg) msg = 'ecliptic longitude, degrees = $' call update(eb, msg) msg = 'year, AD$' call update(year, msg) call ecleq(el, eb, year, alp, del) type *,' right ascension = ',alp,' hours' type *,' declination = ',del,' degrees' go to 10 1100 continue ! here for loc -> solar msg = 'zenith angle, deg = $' call update(theta, msg) msg = 'azimuth (N of E), degrees = $' call update(phi,msg) type *,' enter local time:' msg = 'time zone, +E = $' call update(zone, msg) msg = ' local year = $' call update(year, msg) msg = ' local month = $' call update(month, msg) msg = ' local day = $' call update(day, msg) msg = ' local hour = $' call update(hour, msg) msg = ' local minutes = $' call update(min, msg) msg = ' local seconds = $' call update(sec, msg) msg = ' det lat, deg N = $' call update(detlat, msg) msg = ' det lat, deg W = $' call update(detlong, msg) * get local siderial time call ltlst(year, month, day, hour, min, sec, & detlat, detlong, zone, dlst) msg = 'local siderial time, hrs = $' call update(dlst, msg) * get equatorial coordinates, ra and dec call loceq(theta, phi, dlst, detlat, alp, del) msg = 'got right ascension, hours = $' call update(alp, msg) msg = 'got declination, degrees = $' call update(del, msg) * and then get ecliptic coordinates call eqecl(dalp, del, yr, el, eb) msg = ' got ecliptic longitude, degrees = $' call update(el, msg) msg = ' got ecliptic latitude, degrees = $' call update(eb, msg) * convert local to GMT and dates call ltgmt(dlt, year, month, day, zone, & dgmt, gyear, gmonth, gday) * and get the Julian data at 0 hrs call dayjul(gday, gmonth, gyear, rjd) msg = 'got GMT = $' call update(dgmt, msg) msg = 'got Jul day at 0 hrs G = $' call update(rjd, msg) * now finally get sun's ecliptic coordinates call sunecl(rjd, dgmt, el_sun, sdist) type *,' solar ecliptic latitude = 0.0 degrees (always)' msg = ' solar ecliptic longitude, degrees = $' call update(el_sun, msg) msg = ' earth-sun dist, lt hours = $' call update(sdist, msg) sb = eb ! latitude same for ecliptic and solar coords sl = el - el_sun ! get longitude re sun if(sl .gt. 360.0) sl = sl - 360.0 if(sl .lt. 0.0) sl = sl + 360.0 type *,' long re sun = ',sl,' degrees' * and get the angular distance from point to sun delsl = sl/15.0 ! in hours sb_sun = 0.0 call angdist(delsl, sb, sb_sun, eta) type *,' and dist from sun, deg =',eta go to 10 1200 continue ! here for solar -> loc msg = 'solar longitude, degrees = $' C call update(sl, msg) C msg = 'solar lattitude, degrees = $' call update(sb, msg) C msg = 'decimal local time, hrs = $' call update(dlt, msg) C call solloc(sl, sb, dlt, theta, phi) type *,' sorry not available yet' C type *,' zenith angle = ',theta,' degrees' C type *,' azimuth (N of E) = ',phi,' degrees' go to 10 1300 continue ! here for getting lunar coords msg = ' jd at 0 hrs gmt = $' call update(rjd, msg) msg = ' dgmt, hrs = $' call update(dgmt, msg) call lunar(rjd, dgmt, alp_moon, del_moon) msg = ' got lunar ra, hrs = $' call update(alp_moon, msg) msg = ' got lunar dec, deg = $' call update(del_moon, msg) go to 10 1400 continue ! here for getting angle from moon msg = ' event ra, hrs = $' call update(alp, msg) msg = ' event dec, deg = $' call update(del, msg) msg = ' lunar ra, hrs = $' call update(alp_moon, msg) msg = ' lunar dec, deg = $' call update(del_moon, msg) delalp = alp - alp_moon call angdist(delalp, del, del_moon, eta) msg = ' got angular separation, deg = $' call update(eta, msg) go to 10 1500 continue * set or change nominal source and detector coordinates msg = 'det lat, deg north$' call update(detlat, msg) msg = 'det long, deg west$' call update(detlong, msg) type *,' right ascension of source (hr,min,sec)' call dhhms(alp_src, shr, smin, ssec) msg = ' ra hours$' call update(shr, msg) msg = ' ra minutes$' call update(smin, msg) msg = ' ra seconds$' call update(ssec,msg) call hmsdh(shr, smin, ssec, alp_src) msg = ' ra of source in hours $' call update(alp_src, msg) type *,' declination of source (deg,min,sec)' call dhhms(del_src, ddeg, dmin, dsec) msg = 'dec degrees $' call update(ddeg,msg) msg = 'dec minutes $' call update(dmin,msg) msg = 'dec seconds $' call update(dsec,msg) * note slightly funny call here, but OK (degrees, not hours) call hmsdh(ddeg,dmin,dsec,del_src) msg = ' dec of source in degrees $' call update (del_src,msg) type *,' enter detector local time' msg = 'time zone, +E $' call update(zone,msg) msg = ' local year$' call update(year,msg) msg = ' local month$' call update(month,msg) msg = ' local day$' call update(day,msg) msg = ' local hour$' call update(hour,msg) msg = ' local minutes$' call update(min,msg) msg = ' local seconds$' call update(sec,msg) call hmsdh(hour,min,sec,dlt) msg = ' got decimal local time, hrs $' call update(dlt,msg) call ltgmt(dlt, year, month, day, zone, & dgmt, gyear, gmonth, gday) msg = 'got greenwich year$' call update(gyear, msg) msg = 'got greenwich month$' call update(gmonth, msg) msg = 'got greenwich day$' call update(gday, msg) msg = 'got greenwich hour, dgmt $' call update(dgmt, msg) call dayjul(gday, gmonth, gyear, rjd) msg = ' Julian date at 0 hours Greenwich $' call update(rjd, msg) jd = rjd + dgmt/24.0 msg = ' decimal Julian date Greenwich $' call update(jd, msg) go to 10 1600 continue !barycentric correction msg = ' ra, hours $' call update(alp, msg) msg = ' dec, degrees $' call update (del, msg) msg = 'time zone, +E $' call update(zone, msg) msg = ' local time, hrs $' call update(dlt, msg) msg = ' local year = $' call update(year, msg) msg = ' local month = $' call update(month, msg) msg = ' local day = $' call update(day, msg) type *,' convering to Greenwich time' * get the Greenwich time and date call ltgmt(dlt, year, month, day, zone, & dgmt, gyear, gmonth, gday) msg = ' got dgmt, hrs $' call update(dgmt, msg) type *,' getting Julian day' * get the Julian date of beginning of G day at 0 hrs, rjd call dayjul(gday, gmonth, gyear, rjd) msg = ' got jul day at zero hrs GMT $' call update(rjd, msg) type *,' getting barycenter time corr and sun coords' * and then get the barycenter time correction and sun RA and DEC call barycen(rjd, dgmt, alp, del, & deltat, alp_sun, del_sun, el_sun) type *,' got RA of sun = ',alp_sun,' hours' type *,' got DEC sun = ',del_sun,' degrees' type *,' time correction = ',deltat,' in hours' deltatm = deltat * 60.0 type *,' same correction = ',deltatm,' in min' deltats = deltat * 3600. type *,' same correction = ',deltats,' in sec' type *,' got sun longitude = ',el_sun,' deg along ecliptic' goto 10 * here to convert event coords from mine to ra and dec * and to get angular dist from source 1700 continue msg = 'x dir cos$' call update(xx,msg) msg = 'y dir cos$' call update(yy,msg) msg = 'z dir cos$' call update(zz,msg) fnorm=sqrt(xx**2+yy**2+zz**2) type *, ' norm =',fnorm xx=xx/fnorm yy=yy/fnorm zz=zz/fnorm phi=atan2d(yy,xx) theta=acosd(zz) if(phi .lt. 0.0) phi = phi + 360.0 msg = 'zenith angle re detector, deg$' call update(theta, msg) msg = 'azimuth re det, deg n of e$' call update(phi,msg) msg = 'local sid time, decimal hrs$' call update(dlst, msg) call loceq(theta, phi, dlst, detlat, alp, del) type *,' right ascension = ',alp,' hours' type *,' declination = ',del,' degrees' delalp = alp -alp_src call angdist(delalp, del, del_src, eta) type *,' and dist from source, deg =',eta go to 10 end