Loading plio/date/marstime.py 0 → 100644 +273 −0 Original line number Diff line number Diff line ''' Borrowed from Dave Pawlowski. (http://marstiming.readthedocs.io/en/latest/_modules/marstiming.html) Mars timing information based on MARS24: http://www.giss.nasa.gov/tools/mars24/help/algorithm.html Contains several functions for calculating Mars time parameters from Earth time and vice versa. Probably the most useful functions are: getMTfromTime: gets mars time data given a 6 element time list getSZAfromTime: gets the SZA from a 6 element time list and coordinates getLTfromTime: gets the LTST from a 6 element time list and longitude getUTCfromLS: Estimates the Earth time from LS and a Mars year ''' import datetime from numpy import pi, floor,array,shape, cos, sin,ceil,arcsin,arccos,arange,abs from collections import namedtuple d2R = pi/180. def getJD(iTime): '''Get the Julian date in seconds''' offset = 2440587.5 #JD on 1/1/1970 00:00:00 year = iTime[0] month = iTime[1] day = iTime[2] hour = iTime[3] minute = iTime[4] sec = iTime[5] date = datetime.datetime(year,month,day,hour,minute,sec) iTime = [1970,1,1,0,0,0] year = iTime[0] month = iTime[1] day = iTime[2] hour = iTime[3] minute = iTime[4] sec = iTime[5] ref = datetime.datetime(year,month,day,hour,minute,sec) deltaTime = (date-ref) return deltaTime.total_seconds()/86400. + offset def getUTC(jd): '''Get UTC given jd''' offset = 2440587.5 #JD on 1/1/1970 00:00:00 iTime = [1970,1,1,0,0,0] year = iTime[0] month = iTime[1] day = iTime[2] hour = iTime[3] minute = iTime[4] sec = iTime[5] d1970 = datetime.datetime(year,month,day,hour,minute,sec) return d1970 + datetime.timedelta(seconds=((jd-offset)*86400.)) def getJ2000(iTime): '''get date in J2000 epoch.''' jd = getJD(iTime) T = (jd - 2451545.0)/36525 if iTime[0] < 1972 else 0 conversion = 64.184 + 59* T - 51.2* T**2 - 67.1* T**3 - 16.4* T**4 #convert to Terrestrial Time jdTT = jd+(conversion/86400) return jdTT - 2451545.0 def getMarsParams(j2000): '''Mars time parameters''' Coefs = array( [[0.0071,2.2353,49.409], [0.0057,2.7543,168.173], [0.0039,1.1177,191.837], [0.0037,15.7866,21.736], [0.0021,2.1354,15.704], [0.0020,2.4694,95.528], [0.0018,32.8493,49.095]]) dims = shape(Coefs) #Mars mean anomaly: M = 19.3870 + 0.52402075 * j2000 #angle of Fiction Mean Sun alpha = 270.3863 + 0.52403840*j2000 #Perturbers PBS = 0 for i in range(dims[0]): PBS += Coefs[i,0]*cos(((0.985626* j2000 / Coefs[i,1]) + Coefs[i,2])*d2R) #Equation of Center vMinusM = ((10.691 + 3.0e-7 *j2000)*sin(M*d2R) + 0.623*sin(2*M*d2R) + 0.050*sin(3*M*d2R) + 0.005*sin(4*M*d2R) + 0.0005*sin(5*M*d2R) + PBS) return M, alpha, PBS, vMinusM def getMTfromTime(iTime): '''Get Mars time information. :param iTime: 6 element time list [y,m,d,h,m,s] :returns: a named tuple containing the LS value as well as several parameters necessary for other calculations ''' if isinstance(iTime, datetime.datetime): iTime = [iTime.year, iTime.month, iTime.day, iTime.hour, iTime.minute, iTime.second] DPY = 686.9713 refTime = [1955,4,11,10,56,0] #Mars year 1 rDate = getJD(refTime) thisTime = getJD(iTime) year = floor((thisTime - rDate)/DPY)+1 j2000 = getJ2000(iTime) M,alpha,PBS,vMinusM = getMarsParams(j2000) LS = (alpha + vMinusM) while LS > 360: LS -= 360 if LS < 0: LS = 360. + 360.*(LS/360. - ceil(LS/360.0)) EOT = 2.861*sin(2*LS*d2R)-0.071*sin(4*LS*d2R)+0.002*sin(6*LS*d2R)-vMinusM MTC = (24*(((j2000 - 4.5)/1.027491252)+44796.0 - 0.00096 )) % 24 subSolarLon = ((MTC+EOT*24/360.)*(360/24.)+180) % 360 solarDec = (arcsin(0.42565*sin(LS*d2R))/d2R+0.25*sin(LS*d2R)) data = namedtuple('data','ls year M alpha PBS vMinusM MTC EOT subSolarLon solarDec') d1 = data(ls = LS,year=year,M=M,alpha=alpha,PBS=PBS,vMinusM=vMinusM,MTC=MTC,EOT=EOT, subSolarLon=subSolarLon,solarDec=solarDec) return d1 def getUTCfromLS(marsyear,LS): '''Get a UTC starting with an estimate of LS using an orbit angle approximation then iteratively closing in on the correct LS by incrementing the a day first and then hour. :param marsyear: an int mars year :param ls: ls- mars solar longitude :returns: UTC1 (python datetime)''' #Get LS to within this value: error = 0.001 DPY = 686.9713 ###Start with estimate refTime = [1955,4,11,10,56,0] #Mars year 1 rDate = getJD(refTime) iTime = getUTC(rDate+(marsyear-1)*DPY) #LS 0 of given mars year #Now we have a guess, iterate over the day to get closer and closer. thisTime = [iTime.year,iTime.month,iTime.day,iTime.hour,iTime.minute,iTime.second] thisLS = 0 factor = 1 #do we increment up or down? iTry = 0 dt = 60 #hours. This will get smaller as we get closer counter = 0 olddiff = 1000. diff = 100 while diff > error: iTime = iTime+factor*datetime.timedelta(hours=dt) thisTime = [iTime.year,iTime.month,iTime.day,iTime.hour,iTime.minute,iTime.second] timedata = getMTfromTime(thisTime) thisLS,myear = timedata.ls, timedata.year diff = abs(thisLS - LS) if diff > olddiff: factor = -1*factor counter += 1 if counter > 1: dt = dt/60. olddiff = diff iTry += 1 if iTry > 1000: print('Problem getting UTC from Ls in 2nd diff loop') print('Quitting if function getUTCfromLS...') exit(1) return iTime def getSZAfromTime(time,lon,lat): '''Get SZA from Earth time and Mars coordinates. :param the time: [y,m,d,h,m,s] or datetime object :param lon: the longitude in degrees :param lat: the latitude in degrees :returns: the solar zenith angle (float)''' if type(time) == datetime.datetime: time = [time.year,time.month,time.day,time.hour,time.minute,time.second] timedata = getMTfromTime(time) SZA = arccos(sin(timedata.solarDec*d2R)*sin(lat*d2R)+ cos(timedata.solarDec*d2R)*cos(lat*d2R)*cos((lon-timedata.subSolarLon)*d2R))/d2R return SZA def SZAGetTime(sza,date, lon, lat): '''Find the time on a given date and location when the SZA is a given value. :param sza: Solar zenith angle in degrees :param date: [y,m,d]< :param lon: the longitude in degrees :param lat: the latitude in degrees :returns: A python datetime object ''' thisDate = datetime.datetime(date[0],date[1],date[2]) count = 0 counter = 0 error = 1 factor = 1 dt = 15 #minutes thisSza = getSZAfromTime(thisDate,lon,lat) diff = abs(thisSza - sza) while diff > error: thisDate += factor*datetime.timedelta(minutes=dt) thisSza = getSZAfromTime(thisDate,lon,lat) newdiff = abs(thisSza - sza) if newdiff > diff: factor = -1*factor counter += 1 if counter > 1: #Wait until counter is > 1 in case we start off going the wrong way! dt = dt/2. count += 1 if abs(diff - newdiff)/2. < error and counter > 5: print('this location doesnt reach the given SZA. Returning closest value... {:f}'.format(thisSza)) return thisDate, thisSza diff = newdiff return thisDate, thisSza def getLTfromTime(iTime,lon): '''The mars local solar time from an earth time and mars longitude. :param iTime: 6 element list: [y,m,d,h,m,s] :param lon: the longitude in degrees :returns: The local time (float)''' timedata = getMTfromTime(iTime) LMST = timedata.MTC-lon*(24/360.) LTST = LMST + timedata.EOT*(24/360.) return LTST plio/date/tests/test_marstime.py 0 → 100644 +61 −0 Original line number Diff line number Diff line import numpy as np import unittest import datetime from .. import marstime class TestMarsTime(unittest.TestCase): def equality_test(self): earth = marstime.getUTCfromLS(24, 130) mars = marstime.getMTfromTime(earth) assert(earth.date() == marstime.getUTCfromLS(mars.year,mars.ls).date()) def test_J2000(self): iTime = [2001,11,13,2,45,2] testJD = marstime.getJ2000(iTime) callibration = 58891502.000000 #test should be this value. diff = testJD - callibration print(testJD) print(diff) self.assertTrue(diff == -58890820.38465065) def test_UTC(self): jd = 2452226.614606 date = marstime.getUTC(jd) callibration = datetime.datetime(2001,11,13,2,45,2) diff = callibration - date self.assertTrue(str(diff) == '0:00:00.041599') def test_LS(self): iTime = [2000,1,6,0,0,0] lsdata = marstime.getMTfromTime(iTime) ls = lsdata.ls year = lsdata.year callibration = 277.18677 diff = ls - callibration print(diff) self.assertAlmostEqual(diff, 0.0) def test_SZA(self): '''test getSZAfromTime''' itime = [2000,1,6,0,0,0] lon = 0.0 lat = 0.0 expected = 154.261908004 sza = marstime.getSZAfromTime(itime,lon,lat) self.assertAlmostEqual(sza, expected) def test_SZAGetTime(self): out = marstime.SZAGetTime(0, [1,1,1], 0,0) self.assertAlmostEqual(out[1], 9.454213735250395) def test_LTfromTime(self): '''test getLTfromTime function''' iTime = [2000,1,6,0,0,0] lon = 0.0 LTST = marstime.getLTfromTime(iTime,lon) expected = 23.648468934 self.assertAlmostEqual(LTST, expected) Loading
plio/date/marstime.py 0 → 100644 +273 −0 Original line number Diff line number Diff line ''' Borrowed from Dave Pawlowski. (http://marstiming.readthedocs.io/en/latest/_modules/marstiming.html) Mars timing information based on MARS24: http://www.giss.nasa.gov/tools/mars24/help/algorithm.html Contains several functions for calculating Mars time parameters from Earth time and vice versa. Probably the most useful functions are: getMTfromTime: gets mars time data given a 6 element time list getSZAfromTime: gets the SZA from a 6 element time list and coordinates getLTfromTime: gets the LTST from a 6 element time list and longitude getUTCfromLS: Estimates the Earth time from LS and a Mars year ''' import datetime from numpy import pi, floor,array,shape, cos, sin,ceil,arcsin,arccos,arange,abs from collections import namedtuple d2R = pi/180. def getJD(iTime): '''Get the Julian date in seconds''' offset = 2440587.5 #JD on 1/1/1970 00:00:00 year = iTime[0] month = iTime[1] day = iTime[2] hour = iTime[3] minute = iTime[4] sec = iTime[5] date = datetime.datetime(year,month,day,hour,minute,sec) iTime = [1970,1,1,0,0,0] year = iTime[0] month = iTime[1] day = iTime[2] hour = iTime[3] minute = iTime[4] sec = iTime[5] ref = datetime.datetime(year,month,day,hour,minute,sec) deltaTime = (date-ref) return deltaTime.total_seconds()/86400. + offset def getUTC(jd): '''Get UTC given jd''' offset = 2440587.5 #JD on 1/1/1970 00:00:00 iTime = [1970,1,1,0,0,0] year = iTime[0] month = iTime[1] day = iTime[2] hour = iTime[3] minute = iTime[4] sec = iTime[5] d1970 = datetime.datetime(year,month,day,hour,minute,sec) return d1970 + datetime.timedelta(seconds=((jd-offset)*86400.)) def getJ2000(iTime): '''get date in J2000 epoch.''' jd = getJD(iTime) T = (jd - 2451545.0)/36525 if iTime[0] < 1972 else 0 conversion = 64.184 + 59* T - 51.2* T**2 - 67.1* T**3 - 16.4* T**4 #convert to Terrestrial Time jdTT = jd+(conversion/86400) return jdTT - 2451545.0 def getMarsParams(j2000): '''Mars time parameters''' Coefs = array( [[0.0071,2.2353,49.409], [0.0057,2.7543,168.173], [0.0039,1.1177,191.837], [0.0037,15.7866,21.736], [0.0021,2.1354,15.704], [0.0020,2.4694,95.528], [0.0018,32.8493,49.095]]) dims = shape(Coefs) #Mars mean anomaly: M = 19.3870 + 0.52402075 * j2000 #angle of Fiction Mean Sun alpha = 270.3863 + 0.52403840*j2000 #Perturbers PBS = 0 for i in range(dims[0]): PBS += Coefs[i,0]*cos(((0.985626* j2000 / Coefs[i,1]) + Coefs[i,2])*d2R) #Equation of Center vMinusM = ((10.691 + 3.0e-7 *j2000)*sin(M*d2R) + 0.623*sin(2*M*d2R) + 0.050*sin(3*M*d2R) + 0.005*sin(4*M*d2R) + 0.0005*sin(5*M*d2R) + PBS) return M, alpha, PBS, vMinusM def getMTfromTime(iTime): '''Get Mars time information. :param iTime: 6 element time list [y,m,d,h,m,s] :returns: a named tuple containing the LS value as well as several parameters necessary for other calculations ''' if isinstance(iTime, datetime.datetime): iTime = [iTime.year, iTime.month, iTime.day, iTime.hour, iTime.minute, iTime.second] DPY = 686.9713 refTime = [1955,4,11,10,56,0] #Mars year 1 rDate = getJD(refTime) thisTime = getJD(iTime) year = floor((thisTime - rDate)/DPY)+1 j2000 = getJ2000(iTime) M,alpha,PBS,vMinusM = getMarsParams(j2000) LS = (alpha + vMinusM) while LS > 360: LS -= 360 if LS < 0: LS = 360. + 360.*(LS/360. - ceil(LS/360.0)) EOT = 2.861*sin(2*LS*d2R)-0.071*sin(4*LS*d2R)+0.002*sin(6*LS*d2R)-vMinusM MTC = (24*(((j2000 - 4.5)/1.027491252)+44796.0 - 0.00096 )) % 24 subSolarLon = ((MTC+EOT*24/360.)*(360/24.)+180) % 360 solarDec = (arcsin(0.42565*sin(LS*d2R))/d2R+0.25*sin(LS*d2R)) data = namedtuple('data','ls year M alpha PBS vMinusM MTC EOT subSolarLon solarDec') d1 = data(ls = LS,year=year,M=M,alpha=alpha,PBS=PBS,vMinusM=vMinusM,MTC=MTC,EOT=EOT, subSolarLon=subSolarLon,solarDec=solarDec) return d1 def getUTCfromLS(marsyear,LS): '''Get a UTC starting with an estimate of LS using an orbit angle approximation then iteratively closing in on the correct LS by incrementing the a day first and then hour. :param marsyear: an int mars year :param ls: ls- mars solar longitude :returns: UTC1 (python datetime)''' #Get LS to within this value: error = 0.001 DPY = 686.9713 ###Start with estimate refTime = [1955,4,11,10,56,0] #Mars year 1 rDate = getJD(refTime) iTime = getUTC(rDate+(marsyear-1)*DPY) #LS 0 of given mars year #Now we have a guess, iterate over the day to get closer and closer. thisTime = [iTime.year,iTime.month,iTime.day,iTime.hour,iTime.minute,iTime.second] thisLS = 0 factor = 1 #do we increment up or down? iTry = 0 dt = 60 #hours. This will get smaller as we get closer counter = 0 olddiff = 1000. diff = 100 while diff > error: iTime = iTime+factor*datetime.timedelta(hours=dt) thisTime = [iTime.year,iTime.month,iTime.day,iTime.hour,iTime.minute,iTime.second] timedata = getMTfromTime(thisTime) thisLS,myear = timedata.ls, timedata.year diff = abs(thisLS - LS) if diff > olddiff: factor = -1*factor counter += 1 if counter > 1: dt = dt/60. olddiff = diff iTry += 1 if iTry > 1000: print('Problem getting UTC from Ls in 2nd diff loop') print('Quitting if function getUTCfromLS...') exit(1) return iTime def getSZAfromTime(time,lon,lat): '''Get SZA from Earth time and Mars coordinates. :param the time: [y,m,d,h,m,s] or datetime object :param lon: the longitude in degrees :param lat: the latitude in degrees :returns: the solar zenith angle (float)''' if type(time) == datetime.datetime: time = [time.year,time.month,time.day,time.hour,time.minute,time.second] timedata = getMTfromTime(time) SZA = arccos(sin(timedata.solarDec*d2R)*sin(lat*d2R)+ cos(timedata.solarDec*d2R)*cos(lat*d2R)*cos((lon-timedata.subSolarLon)*d2R))/d2R return SZA def SZAGetTime(sza,date, lon, lat): '''Find the time on a given date and location when the SZA is a given value. :param sza: Solar zenith angle in degrees :param date: [y,m,d]< :param lon: the longitude in degrees :param lat: the latitude in degrees :returns: A python datetime object ''' thisDate = datetime.datetime(date[0],date[1],date[2]) count = 0 counter = 0 error = 1 factor = 1 dt = 15 #minutes thisSza = getSZAfromTime(thisDate,lon,lat) diff = abs(thisSza - sza) while diff > error: thisDate += factor*datetime.timedelta(minutes=dt) thisSza = getSZAfromTime(thisDate,lon,lat) newdiff = abs(thisSza - sza) if newdiff > diff: factor = -1*factor counter += 1 if counter > 1: #Wait until counter is > 1 in case we start off going the wrong way! dt = dt/2. count += 1 if abs(diff - newdiff)/2. < error and counter > 5: print('this location doesnt reach the given SZA. Returning closest value... {:f}'.format(thisSza)) return thisDate, thisSza diff = newdiff return thisDate, thisSza def getLTfromTime(iTime,lon): '''The mars local solar time from an earth time and mars longitude. :param iTime: 6 element list: [y,m,d,h,m,s] :param lon: the longitude in degrees :returns: The local time (float)''' timedata = getMTfromTime(iTime) LMST = timedata.MTC-lon*(24/360.) LTST = LMST + timedata.EOT*(24/360.) return LTST
plio/date/tests/test_marstime.py 0 → 100644 +61 −0 Original line number Diff line number Diff line import numpy as np import unittest import datetime from .. import marstime class TestMarsTime(unittest.TestCase): def equality_test(self): earth = marstime.getUTCfromLS(24, 130) mars = marstime.getMTfromTime(earth) assert(earth.date() == marstime.getUTCfromLS(mars.year,mars.ls).date()) def test_J2000(self): iTime = [2001,11,13,2,45,2] testJD = marstime.getJ2000(iTime) callibration = 58891502.000000 #test should be this value. diff = testJD - callibration print(testJD) print(diff) self.assertTrue(diff == -58890820.38465065) def test_UTC(self): jd = 2452226.614606 date = marstime.getUTC(jd) callibration = datetime.datetime(2001,11,13,2,45,2) diff = callibration - date self.assertTrue(str(diff) == '0:00:00.041599') def test_LS(self): iTime = [2000,1,6,0,0,0] lsdata = marstime.getMTfromTime(iTime) ls = lsdata.ls year = lsdata.year callibration = 277.18677 diff = ls - callibration print(diff) self.assertAlmostEqual(diff, 0.0) def test_SZA(self): '''test getSZAfromTime''' itime = [2000,1,6,0,0,0] lon = 0.0 lat = 0.0 expected = 154.261908004 sza = marstime.getSZAfromTime(itime,lon,lat) self.assertAlmostEqual(sza, expected) def test_SZAGetTime(self): out = marstime.SZAGetTime(0, [1,1,1], 0,0) self.assertAlmostEqual(out[1], 9.454213735250395) def test_LTfromTime(self): '''test getLTfromTime function''' iTime = [2000,1,6,0,0,0] lon = 0.0 LTST = marstime.getLTfromTime(iTime,lon) expected = 23.648468934 self.assertAlmostEqual(LTST, expected)
