C+ C NAME: C Time2Carrington C PURPOSE: C Calculate the 'modified Carrington variable' for a given spacecraft C position at a given time tt and v.v. C CATEGORY: C Celestial mechanics C CALLING SEQUENCE: subroutine Time2Carrington(id,tt,nc,xc) C INPUTS: C id integer id=0: time --> Carrington variable C id=1: Carington variable --> time C tt(2) integer id=0: UT time C nc integer id=1: Carrington rotation C xc double precision id=1: fraction of Carrington rotation C OUTPUTS: C nc integer id=0: Carrington rotation C xc double precision id=0: fraction of Carrington rotation C tt(2) integer id=1: UT time C INCLUDE: include 'sun.h' include 'fortime.h' C CALLS: C Time2JD, Time2ChangeUnits, Time2ConvertUnits, Time2Delta C HLngEarth, Time2Day8, TimeAdd, Time2to1, Time2Add, Say C RESTRICTIONS: C > It is assumed that a time origin has been defined (using e.g. TimeInit()). C > Calculations are done internally with time in milli-second precision. C This has several side-effects: C - If the time units in effect are smaller then the origin should be put at an integer C multiple of milli-second, or it will be rounded to the nearest milli-second. C - For id=0, the actual time is rounded to the nearest milli-second C - For id=1, the returned time will be rounded to the nearest time unit in effect when C the subroutine was called. C PROCEDURE: C > DEFINITION: Modified Carrington variable C For a given time 'tt' the spacecraft position is inside Carrington C rotation 'nc' at heliographic location (Lng,Lat). C The rotation number 'nc' and the longitude Lng are combined in the C 'modified Carrington variable' by C xc = nc+(360-Lng)/360 C i.e. the integer part of 'xc' identifies the Carrington rotation, C and the fraction the heliographic longitude; as the spacecraft moves C across Carrington rotation 'nc', 'xc' increases from 'nc' to 'nc'+1. C C > The approximate start time is estimated based on the starttime of C Carrington rotation 1750 and the synodic rotation rate of Earth. C A refinement is made using the actual heliographic longitude of Earth C returned by HLngEarth. C MODIFICATION HISTORY: C JAN-2004, Paul Hick (UCSD/CASS) C FEB-2004, Paul Hick (UCSD/CASS) C Fixed bug in iteration loop for conversion from Carrington variable C to time. The iteration would sometimes go into an infinite loop, C bouncing between two times 1 msec apart. C APR-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu) C Fixed bug that disables abort if iteration count is exceeded. C- integer tt(2) integer nc double precision xc character cSay*15 /'Time2Carrington'/ ! Start of Carrington rotation 1750 is JD 2445871.91 ! Fraction of day is stored in msec. integer units /TIME__MSEC_IN_DAY/ integer t1750(2) /2445871,78624000/ integer old_units integer t2 (2) integer dt2(2) integer dt1 integer dt1old double precision dd double precision HLng double precision HLngEarth if (id .eq. 0) then !------- ! xc MUST be based on the output of HLngEarth for tt, not the time ! rounded to nearest milli-second (t2 below). This ensures that ! xc is consistent with the output of HLngEarth very close to the ! edge of a rotation. xc = (360.0d0-HLngEarth(tt))/360.0d0 ! Exact positive fraction !------- ! This changes the origin if units used are smaller than milli-second ! and the origin is defined with sub-millisecond precision. old_units = TIME__ORIGIN(3) call Time2ChangeUnits(units) ! Here round-off to the nearest millisecond takes place. call Time2ConvertUnits(old_units,units,tt,t2) call Time2JD(0,t2,t2) ! Time to Julian days call Time2Delta(t2,t1750,t2) call Time2Day8(0,t2,dd) dd = dble(dd)/SUN__SYNODICP nc = int(dd) dd = dd-dble(nc) ! Approximate positive fraction nc = 1750+nc dd = dd-xc ! Error in fraction [-1,+1] if (dd .lt. -0.5d0) then nc = nc-1 else if (dd .gt. 0.5d0) then nc = nc+1 end if else !------- ! This changes the origin if units used are smaller than milli-second ! and the origin is defined with sub-millisecond precision. old_units = TIME__ORIGIN(3) call Time2ChangeUnits(units) !------- ! Calculate approximate start time of rotation call Time2Day8(1,t2,(xc+dble(nc-1750))*SUN__SYNODICP) call Time2Add(t1750,t2,t2) ! Approx time in Julian days call Time2JD(1,t2,t2) ! Approx time relative to current time origin !------- ! The first difference dd will lie inside [-360,360]. Usually dd will be very ! small (order of a degree or so), but occasionally it may be close to +/- 360 ! (if HLngEarth and HLng are close to the edge of a rotation, but at ! opposite sides of the edge. HLng = 360.0d0*(1.0d0-(xc-int(xc))) ! Target heliographic longitude in [0,360] dd = HLngEarth(t2)-HLng ! [-360,+360] if (dd .lt. -180.0d0) then dd = dd+360.0d0 else if (dd .gt. 180.0d0) then dd = dd-360.0d0 end if !------- ! dd is now very small. Convert to a time measure using the synodic rotation ! period of the Sun. Use it iteratively to improve the time estimate call Time2Day8(1,dt2,dd/360.0d0*SUN__SYNODICP) call Time2to1(dt2,dt1) icount = 0 icount_max = 12 do while (dt1 .ne. 0 .and. icount .lt. icount_max) call Time2Add(t2,dt2,t2) dd = HLngEarth(t2)-HLng ! [-360,+360] if (dd .lt. -180.0d0) then dd = dd+360.0d0 else if (dd .gt. 180.0d0) then dd = dd-360.0d0 end if call Time2Day8(1,dt2,dd/360.0d0*SUN__SYNODICP) dt1old = dt1 call Time2to1(dt2,dt1) ! If dt1old = -dt1 the iteration loop will start bouncing between ! two times. If abs(dt1) = 1 we declare success and exit the loop ! by explicitly setting dt1 = 0. if (abs(dt1) .eq. 1 .and. dt1 .eq. -dt1old) dt1 = 0 icount = icount+1 end do ! A typical count should be 3 or 4. Abort if we actually need icount_max ! (set to three times the expected count). if (icount .ge. icount_max) call Say(cSay,'E','Iteration','count exceeded') !------- ! Change to time units in effect when subroutine was called. ! Round-off will occur here when old_units were larger than millisecond. call Time2ConvertUnits(units,old_units,t2,tt) end if !------- ! Restore input time units for origin call Time2ChangeUnits(old_units) return end