C Computation of Greenwich Mean Sidereal Time (GMST) C Greenwich Apparent Sidereal Time (GAST) C Earth Rotation Angle (ERA) C HL ARI Sept. 2008 C HL ARI verbessert 04.04.11 Delta(T) = 67 sec (2011) beruecksichtigt 04.04.11 C HL ARI verbessert 05.03.13 jeweiliges deltat mit Function angebunden 05.03.13 C Genauigkeit: z.B. Berechnung fuer 1960-07-01T07:21:36.572 05.03.13 C mit richtigem deltat = 33 sec 05.03.13 C GMST = 1 59 0.562461 GAST = 1 59 0.380769 ERA = 30 15 30.110952 05.03.13 C mit falschem deltat = 68 sec (fuer 2015) 05.03.13 C GMST = 1 59 0.562465 GAST = 1 59 0.380771 ERA = 30 15 30.110952 05.03.13 C siehe auch Almanac 2013 B8 05.03.13 C can be checked with Astronomical Almanac B13-B20, B21-B24 implicit real*8 (a-h,o-z) double precision iau_GMST06, iau_GST06A, iau_ERA00 data daydec, zero / 2 * 0.d0 / data fracs / 0.d0 / 21.02.11 C C input: (time in UT1) jahr = 2015 Almanac 2015 B12 monat = 7 itag = 8 ihour = 9 imin = 44 isec = 30 fracs = 0.0 d0 26.02.09 C call iau_CAL2JD ( jahr,monat,itag, DJM0, DJM, J ) if(j.ne.0) then 02.10.08 write(6,*) ' FEHLER bei CAL2JD - AUFRUF' 02.10.08 stop 02.10.08 end if 02.10.08 DJ = DJM0 + DJM daydec = (dble(ihour)*3600.d0 C & + dble(imin)*60.d0 + dble(isec))/86400.d0 & + dble(imin)*60.d0 + dble(isec)+ fracs)/86400.d0 26.02.09 C daydec = 1.d0 C write(6,*) DJ, daydec pi = 4.d0 * atan(1.d0) C C +++++++++++++++++ Nur der folgende Teil muss umgestellt werden ++++++++++++++++++++++ C xtta = DJ vor 04.04.11 C xttb = daydec vor 04.04.11 C xuta = xtta vor 04.04.11 C xuta = xtta + 67.d0/86400.d0 vor 04.04.11 xuta = DJ 04.04.11 xutb = daydec 04.04.11 xtta = xuta 04.04.11 C xttb = daydec + 67.d0/86400.d0 04.04.11 C xttb = daydec + 68.d0/86400.d0 28.02.13 fuer APFS 2015 C xttb = daydec + deltat(jahr)/86400.d0 05.03.13 C dt1 = deltat(jahr) 05.03.13 call iau_DAT (jahr,monat,itag,daydec,dAT,ixflag) 19.03.13 write(6,777) jahr,dat 19.03.13 777 format(/' Fuer das Jahr ',I4,' ist dAT :',F8.1,' sec') 05.03.13 if(ixflag.lt.0) then 19.03.13 write(6,776) ixflag 19.03.13 776 format(/' Fuer das Jahr ',I4,' ist ixflag:',I10 ,' STOP') 19.03.13 end if 19.03.13 if(ixflag.gt.0) then 19.03.13 write(6,775) jahr,ixflag 19.03.13 775 format(/' Fuer das Jahr ',I4,' ist ixflag:',I10 ,' WARNING ') 19.03.13 end if 19.03.13 C 19.03.13 xttb = daydec + (dAT + 32.184d0)/86400.d0 19.03.13 C 19.03.13 write(6,774) jahr,dAT+32.184d0 19.03.13 774 format(/' Fuer das Jahr ',I4,' ist deltat:',F8.1,' sec') 19.03.13 C +++++++++++++++++++ Umstellung ENDE ++++++++++++++++++++++++++++++++++++++++++++++++ write(6,*) xtta, ' xtta' write(6,*) xttb, ' xttb' write(6,*) xuta, ' xuta' write(6,*) xutb, ' xutb' GMST = iau_GMST06 ( xuta,xutb,xtta,xttb) * 180.d0/pi/15.d0 write(6,*) GMST,' GMST' call STMISE (GMST,ihour,min,sec) write(6,887) ihour , min, sec, ' GMST [h, m, s]' GAST = iau_GST06A ( xuta,xutb,xtta,xttb) * 180.d0/pi/15.d0 write(6,*) GAST,' GAST' call STMISE (GAST,ihour,min,sec) write(6,887) ihour , min, sec, ' GAST [h, m, s]' C Achtung: ERA immer nur ganze Tage, keine Bruchteile in DJ einfuehren xtta = DJ xttb = zero xuta = xtta xutb = zero C Faktor 15.041.... aus Almanac B8 (gilt auch noch z.B. fuer Almanac 2012) C fuehrt zu Ungenauigkeiten bis zu 2 mas in 24 h 02.10.08 C ==> times.era-interpolation-accuracy.f 02.10.08 ERA = iau_ERA00 ( xuta,xutb) * 180.d0/pi & + daydec * 24.d0 * 15.0410672 d0 if(ERA.gt.360.d0) ERA = ERA - 360.d0 if(ERA.gt.360.d0) ERA = ERA - 360.d0 if(ERA.gt.360.d0) ERA = ERA - 360.d0 call STMISE (ERA ,ihour,min,sec) write(6,887) ihour , min, sec, & ' ERA [deg, arcmin, arcsec]' 887 format(2i8, F10.6, a) stop end subroutine STMISE (arg,ihour,min,sec) implicit real*8 (a-h,o-z) ihour = int(arg) min = int(abs(arg*60.d0) - dble(iabs(ihour))*60.d0) sec = abs(arg)*3600.d0 & -dble(iabs(ihour))*3600.d0-dble(min)*60.d0 return end C include '/work/Tux3/lenhardt/APFS-SBR/IAU-SOFA-20080301/iau_sbr' C include '/work/Tux3/lenhardt/APFS-SBR/IAU-SOFA-20101201/iau_sbr' C include '/work/Tux3/lenhardt/APFS-SBR/IAU-SOFA-20120710/iau_sbr' 05.03.13 include '/work/Tux3/lenhardt/APFS-SBR/IAU-SOFA-20150209/iau_sbr' 16.02.15