OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
astreduc.py
Go to the documentation of this file.
1 '''
2 Created on Nov 24, 2015
3 
4 @author: rhealy (richard.healy@nasa.gov)
5 
6 '''
7 
8 
9 import sys
10 from datetime import datetime as DT
11 from datetime import timedelta as TD
12 import numpy as np
13 
14 MMM={'JAN':1,'FEB':2,'MAR':3,'APR':4,'MAY':5,'JUN':6,'JUL':7,'AUG':8,'SEP':9,'OCT':10,'NOV':11,'DEC':12}
15 
16 edf_delimiter = [2,2,2,9,3,9,9,1,9,9,3,10,10,1,7,7]
17 edf_usecols = [0,1,2,5,8,11,14]
18 edf_dtype = [np.int32,np.int32,np.int32,np.float64,np.float64,np.float64,np.float64]
19 edf_names = ['year','month','day','pmx','pmy','ut1mutc','loda']
20 
21 ls_delimiter = [6,3,3,5,10,10,12]
22 ls_usecols = [0,1,2,6]
23 ls_dtype = [np.int32,(np.str_,3),np.int32,np.float64]
24 ls_names = ['yfo','mco','dfo','taimutco']
25 
26 rad2Deg = 180.0/np.pi
27 deg2Rad = np.pi/180.0
28 as2Deg = 1.0/3600.0
29 as2Rad = as2Deg*deg2Rad
30 """
31 if __name__ == '__main__':
32  import os
33  try:
34  ocvarroot = os.environ['OCVARROOT']
35  except KeyError:
36  print("OCSSW environement variable OCVARROOT not set.")
37  print("Typically it is set to $OCSSW/run/var.")
38  sys.exit()
39  filein = "{}/hico/{}".format(ocvarroot,'nutation.dat')
40 
41  nut = initReduc(filein)
42 
43  print(nut.head(10))
44 
45 class HicoSec2Date():
46  def __init__(self, delta):
47  # From MJM 2012/10/24 Adjusting the time from the epoch (2009-01-01-00:00:00 UTC).
48  # Input is GPS seconds since the epoch. Output is UTC date/time.
49  # adapted to python by RJH @ NASA 11/24/2015
50  mlen = np.array([31,28,31,30,31,30,31,31,30,31,30,31])
51  dinc=delta
52  idinc=np.floor(delta)
53  #Seconds since epoch of the (only) leap second so far
54  ls_2012_06_30 = 86400*((365*3)+(31+29+31+30+31+30))+1
55  # use epoch definition as start
56  self.yyyy=2009
57  self.mm=1
58  self.dd=1
59  self.hour= np.dtype('int32')
60  self.min = np.dtype('int32')
61  self.hour=0
62  self.min=0
63  self.sec=0.e0
64  if idinc > ls_2012_06_30:
65  #date after leap second, so adjust this to new epoch
66  self.yyyy=2012
67  self.mm=7
68  self.dd=1
69  idinc=idinc-ls_2012_06_30
70  dinc=dinc-ls_2012_06_30
71  self.hour=np.int(0)
72  self.min=np.int(0)
73  self.sec=0.e0
74 
75  mm = self.mm
76  mmlen=mlen[mm - 1]
77  jdinc=np.dtype('int32')
78  while idinc >= 86400*mmlen:
79  if mm == 2:
80  if isLeapYear(self.yyyy) == 1:
81  mmlen = mlen[mm-1]+1
82  else:
83  mmlen = mlen[mm-1]
84  else:
85  mmlen = mlen[mm-1]
86  jdinc = 86400*mmlen
87  idinc=idinc-jdinc
88  dinc=dinc - 86400*mmlen
89  if mm != 12:
90  mm=mm+1
91  else:
92  mm=1
93  self.yyyy=self.yyyy+1
94 
95  self.mm = mm
96 # at this point, we have less a full month. Let's figure out how many whole days
97  nd=np.int(idinc/86400)
98  self.dd=self.dd+nd
99  idinc=idinc - 86400*nd
100  dinc=dinc - 86400*nd
101 # Less than 24 hours left
102  self.hour=np.int(idinc/3600)
103  idinc=idinc - self.hour*3600
104  dinc=dinc - self.hour*3600
105 # Less than 1 hour left
106  self.min=np.int(idinc/60)
107  idinc=idinc - self.min*60
108  dinc=dinc - self.min*60
109 # Less than 1 minute left
110  self.sec=dinc"""
111 
112 def get_astrodate(delta_seconds):
113  d_epoch = DT(year=2009, month=1, day=1)
114  delta_t = TD(seconds=delta_seconds)
115  return d_epoch + delta_t
116 
117 def isLeapYear(yyyy):
118  if (yyyy % 4) == 0:
119  if (yyyy % 100) == 0:
120  leap=1
121  else:
122  if (yyyy % 400) == 0:
123  leap=1
124  else:
125  leap=0
126  else:
127  leap=0
128 
129  return leap
130 
131 def increment_seconds_to_date(delta):
132  # From MJM 2012/10/24 Adjusting the time from the epoch (2009-01-01-00:00:00 UTC).
133  # Input is GPS seconds since the epoch. Output is UTC date/time.
134  # adapted to python by RJH @ NASA 11/24/2015
135  mlen = np.array([31,28,31,30,31,30,31,31,30,31,30,31])
136  dinc=delta
137  idinc=int(np.floor(delta))
138  #Seconds since epoch of the (only) leap second so far
139  ls_2012_06_30 = 86400*((365*3)+(31+29+31+30+31+30))+1
140  #else use epoch definition as start
141  yyyy=2009
142  mm=1
143  dd=1
144  if idinc > ls_2012_06_30:
145  #date after leap second, so adjust this to new epoch
146  yyyy=2012
147  mm=7
148  dd=1
149  idinc=idinc-ls_2012_06_30
150  dinc=dinc-ls_2012_06_30
151 
152  mmlen=mlen[mm - 1]
153  jdinc=np.dtype('int32')
154  while idinc >= 86400*mmlen:
155  if mm == 2:
156  if isLeapYear(yyyy) == 1:
157  mmlen = mlen[mm-1]+1
158  else:
159  mmlen = mlen[mm-1]
160  else:
161  mmlen = mlen[mm-1]
162  jdinc = 86400*mmlen
163  dinc=dinc - 86400*mmlen
164  idinc=idinc-jdinc
165  if mm != 12:
166  mm=mm+1
167  else:
168  mm=1
169  yyyy=yyyy+1
170 # at this point, we have less a full month. Let's figure out how many whole days
171  nd=int(idinc/86400)
172  dd=dd+nd
173  idinc=idinc - 86400*nd
174  dinc=dinc - 86400*nd
175 # Less than 24 hours left
176  hour=int(idinc/3600)
177  idinc=idinc - hour*3600
178  dinc=dinc - hour*3600
179 # Less than 1 hour left
180  minute=int(idinc/60)
181  idinc=idinc - min*60
182  dinc=dinc - min*60
183 # Less than 1 minute left
184  sec=dinc
185 
186  return yyyy,mm,dd,hour,minute,sec
187 
188 def initReduc(filename):
189  import pandas as pd
190  Convrt= 0.0001e0/3600.0e0 # 0.0001 # to deg
191  nutation_data = pd.read_csv(filename,skiprows=110,skipfooter=552-216,skipinitialspace=True,sep=' ',
192  names=['a1','a2','a3','a4','a5','A','B','C','D','NUM'], engine='python')
193  nutation_data['A'] = nutation_data['A']*Convrt
194  nutation_data['B'] = nutation_data['B']*Convrt
195  nutation_data['C'] = nutation_data['C']*Convrt
196  nutation_data['D'] = nutation_data['D']*Convrt
197 
198  return nutation_data
199 
200 def getEDFData(filename,year,mon,day):
201  return np.concatenate([ [row[3]*as2Rad,row[4]*as2Rad,row[5],row[6]/1000.] for row in parseEDFfile(filename) if (row[0] == year and row[1] == mon and row[2] == day )])
202 
203 def parseEDFfile(filename,skiphdr=0):
204  return np.genfromtxt(filename,skip_header=skiphdr,
205  delimiter = edf_delimiter,
206  usecols = edf_usecols,
207  dtype = edf_dtype,
208  names = edf_names)
209 
210 def getLSData(filename,year,mon,day):
211  mdate=year*10000 + mon*100 + day
212  lsdat = np.concatenate([ [row[3]] for row in parseLSfile(filename,1) if (mdate > (row[0]*10000 + MMM[row[1]]*100 + row[2] ) )])
213 
214  return lsdat[-1]
215 
216 def parseLSfile(filename,skiphdr=0):
217  return np.genfromtxt(filename,skip_header=skiphdr,
218  delimiter = ls_delimiter,
219  usecols = ls_usecols,
220  dtype = ls_dtype,
221  names = ls_names)
222 
223 ''''
224 converted to python by R. Healy 11/27/2015
225 * ----------------------------------------------------------------------------
226 *
227 * SUBROUTINE POLARM
228 *
229 * this function calulates the transformation matrix for polar motion.
230 * the units for polar motion are input in rad because it's more
231 * efficient to do this in the main routine.
232 *
233 * Author : David Vallado 719-573-2600 1 Mar 2001
234 *
235 * Inputs Description Range / Units
236 * xp - Polar motion coefficient rad
237 * yp - Polar motion coefficient rad
238 *
239 * Outputs :
240 * PM - Polar motion transformation (pef-ecef)
241 *
242 * Locals :
243 * None.
244 *
245 * Coupling :
246 * ROT2 - Rotation about the second axis
247 * ROT3 - Rotation about the third axis
248 *
249 * References :
250 * Vallado 2007, 228
251 *
252 * ----------------------------------------------------------------------------
253 '''
254 def polarm(xp,yp):
255 
256 
257  cosxp = np.cos(xp)
258  sinxp = np.sin(xp)
259  cosyp = np.cos(yp)
260  sinyp = np.sin(yp)
261 
262  pm = np.zeros((3,3),dtype=np.float64)
263 
264  pm[0][0] = cosxp
265  pm[0][1] = sinxp * sinyp
266  pm[0][2] = sinxp * cosyp
267 
268  pm[1][0] = 0.0
269  pm[1][1] = cosyp
270  pm[1][2] = -sinyp
271 
272  pm[2][0] = -sinxp
273  pm[2][1] = cosxp * sinyp
274  pm[2][2] = cosxp * cosyp
275 
276  return pm
277 
278 def quat_to_rot(q):
279  '''
280  #*****************************************************************************80
281  # Converted to Python by R. Healy 11/30/2015
282  # from q_to_r.f90
283  #
284  # Very slightly modified my Marcos Montes.
285  # MJM Note: scalar is q(1).
286  ## ROTATION_QUAT2MAT_3D converts rotation from quaternion to matrix form in 3D.
287  #
288  # Licensing:
289  #
290  # This code is distributed under the GNU LGPL license.
291  #
292  # Modified:
293  #
294  # 27 July 1999
295  #
296  # Author:
297  #
298  # John Burkardt
299  #
300  # Reference:
301  #
302  # James Foley, Andries van Dam, Steven Feiner, John Hughes,
303  # Computer Graphics, Principles and Practice,
304  # Second Edition,
305  # Addison Wesley, 1990.
306  #
307  # Parameters:
308  #
309  # Input, real ( kind = 8 ) Q(4), the quaternion representing the rotation.
310  #
311  # Output, real ( kind = 8 ) a[3,3), the rotation matrix.
312  #
313  '''
314  sin_phi = np.sqrt ( np.sum ( np.square(q[1:4]) ) )
315 
316  cos_phi = q[0]
317 
318  angle = 2.0 * np.arctan2( sin_phi, cos_phi )
319 
320  if sin_phi == 0.0:
321  v1 = 1.0
322  v2 = 0.0
323  v3 = 0.0
324  else:
325  v1 = q[1] / sin_phi
326  v2 = q[2] / sin_phi
327  v3 = q[3] / sin_phi
328 
329  ca = np.cos ( angle )
330  sa = np.sin ( angle )
331 
332  a = np.zeros((3,3))
333  a[0,0] = v1 * v1 + ca * ( 1.0 - v1 * v1 )
334  a[0,1] = ( 1.0 - ca ) * v1 * v2 - sa * v3
335  a[0,2] = ( 1.0 - ca ) * v1 * v3 + sa * v2
336 
337  a[1,0] = ( 1.0 - ca ) * v2 * v1 + sa * v3
338  a[1,1] = v2 * v2 + ca * ( 1.0 - v2 * v2 )
339  a[1,2] = ( 1.0 - ca ) * v2 * v3 - sa * v1
340 
341  a[2,0] = ( 1.0 - ca ) * v3 * v1 - sa * v2
342  a[2,1] = ( 1.0 - ca ) * v3 * v2 + sa * v1
343  a[2,2] = v3 * v3 + ca * ( 1.0 - v3 * v3 )
344 
345  return a
346 
347 
348 def hms2ut(Hr, Minute, Sec):
349  '''
350  * Converted to python by R. Healy 11/30/2015
351  * broken into two functions from fortran sub HMS_SEC:
352  hms2ut
353  ut2hms
354  * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
355  *
356  * SUBROUTINE HMS_SEC
357  *
358  * this subroutine converts Hours, Minutes and Seconds into seconds from the
359  * beginning of the day.
360  *
361  * Author : David Vallado 719 - 573 - 2600 1 Mar 2001
362  *
363  * Inputs Description Range / Units
364  * Hr - Hours 0 .. 24
365  * minute - Minutes 0 .. 59
366  * Sec - Seconds 0.0D0 .. 59.99D0
367  * Direction - Which set of vars to output FROM TOO
368  *
369  * OutPuts :
370  * Sec - Seconds 0.0D0 .. 86400.0D0
371  *
372  * Locals :
373  * Temp - Temporary variable
374  *
375  * Coupling :
376  * None.
377  *
378  * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
379 '''
380  UTSec= Hr*3600.0 + Minute*60.0 + Sec
381 
382  return UTSec
383 
384 def ut2hms(UTSec):
385  '''
386  * Converted to python by R. Healy 11/30/2015
387  * broken into two functions from fortran sub HMS_SEC:
388  hms2ut
389  ut2hms
390  * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
391  *
392  * SUBROUTINE HMS_SEC
393  *
394  * this subroutine converts Hours, Minutes and Seconds into seconds from the
395  * beginning of the day.
396  *
397  * Author : David Vallado 719 - 573 - 2600 1 Mar 2001
398  *
399  * Inputs Description Range / Units
400  * Hr - Hours 0 .. 24
401  * minute - Minutes 0 .. 59
402  * Sec - Seconds 0.0D0 .. 59.99D0
403  * Direction - Which set of vars to output FROM TOO
404  *
405  * OutPuts :
406  * Sec - Seconds 0.0D0 .. 86400.0D0
407  *
408  * Locals :
409  * Temp - Temporary variable
410  *
411  * Coupling :
412  * None.
413  *
414  * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
415 '''
416 
417  Temp= UTSec / 3600.0
418  Hr = np.floor( Temp )
419  Minute = np.floor( (Temp - Hr)*60.0 )
420  Sec = (Temp - Hr - Minute/60.0 ) * 3600.0
421 
422  return Hr, Minute, Sec
423 
424 
425 def JDay ( Year,Mon,Day,Hr,minute, Sec ):
426  '''
427  * Converted to python by R. Healy
428  * -----------------------------------------------------------------------------
429  *
430  * SUBROUTINE JDay
431  *
432  * this subroutine finds the Julian date given the Year, Month, Day, and Time.
433  *
434  * Author : David Vallado 719-573-2600 1 Mar 2001
435  *
436  * Inputs Description Range / Units
437  * Year - Year 1900 .. 2100
438  * Mon - Month 1 .. 12
439  * Day - Day 1 .. 28,29,30,31
440  * Hr - Universal Time Hour 0 .. 23
441  * minute - Universal Time minute 0 .. 59
442  * Sec - Universal Time Sec 0.0D0 .. 59.999D0
443  * WhichType - Julian .or. Gregorian calender 'J' .or. 'G'
444  *
445  * Outputs :
446  * JD - Julian Date days from 4713 BC
447  *
448  * Locals :
449  * B - Var to aid Gregorian dates
450  *
451  * Coupling :
452  * None.
453  *
454  * References :
455  * Vallado 2007, 189, Alg 14, Ex 3-14
456  * -----------------------------------------------------------------------------
457  '''
458 
459  jday = 367.0 * Year \
460  - np.int( (7* (Year+np.int ( (Mon+9)/12) ) ) * 0.25 ) \
461  + np.int( 275*Mon / 9 ) \
462  + Day + 1721013.5 \
463  + ( (Sec/60.0 + minute ) / 60.0 + Hr ) / 24.0
464  return jday
465 
466 class UT2time():
467  '''
468  * Converted to python by R. Healy 11/30/2015
469  * ------------------------------------------------------------------------------
470  *
471  * SUBROUTINE CONVTIME
472  *
473  * this subroutine finds the time parameters and Julian century values for inputs
474  * of UTC or UT1. Numerous outputs are found as shown in the local variables.
475  * Because calucations are in UTC, you must include TimeZone IF ( you enter a
476  * local time, otherwise it should be zero.
477  *
478  * Algorithm : A file of record contains the timing data
479  * Seeks are performed to obtain the data
480  * Data starts Jan 1, 1980, thus JD = 2444238.5D0 in the code
481  * Calculate the answer depending on initial time type
482  *
483  * Author : David Vallado 719-573-2600 1 Mar 2001
484  *
485  * Inputs Description Range / Units
486  * Year - Year 1900 .. 2100
487  * Mon - Month 1 .. 12
488  * Day - Day 1 .. 28,29,30,31
489  * Hr - Universal Time Hour 0 .. 23
490  * minute - Universal Time minute 0 .. 59
491  * SEC - Universal Time SEC 0.0D0 .. 59.999D0
492  * TimeZone - Offset to UTC from local SITE 0 .. 23 hr
493  * TypeUTIn - Type of input UT 1 (UT1), else UTC
494  * DUT1 - Delta of UTC - UT1 SEC
495  *
496  * Outputs :
497  * DAT - Delta of TAI-UTC SEC [MJM: THIS IS AN INPUT]
498  * xp - Polar motion coefficient arcsec [MJM: INPUT]
499  * yp - Polar motion coefficient arcsec [MJM: INPUT]
500  * UT1 - Universal time SEC
501  * TUT1 - Julian centuries of UT1
502  * JDUT1 - Julian Date of UT1 days from 4713 BC
503  * UTC - Coordinated Universal Time SEC
504  * TAI - Atomic time SEC
505  * TDT - Terrestrial Dynamical time SEC
506  * TTDT - Julian centuries of TDT
507  * JDTDT - Julian Date of TDT days from 4713 BC
508  * TDB - Terrestrial Barycentric time SEC
509  * TTDB - Julian centuries of TDB
510  * JDTDB - Julian Date of TDB days from 4713 BC
511  *
512  * Locals :
513  * HrTemp - Temporary hours hr
514  * MinTemp - Temporary miNutes minute
515  * SecTemp - Temporary seconds SEC
516  * LocalHr - Difference to local time hr
517  * JD - Julian Date of request days from 4713 BC
518  * ME - Mean Anomaly of the Earth rad
519  * TimeFile - File of record with time data
520  * CurrTimeRec - Current Time record
521  *
522  * Coupling :
523  * HMS_SEC - Conversion between hr-minute-SEC .and. seconds
524  * jday - Find the Julian date
525  *
526  * References :
527  * vallado 2007, 201, alg 16, ex 3-7
528  *
529  * ------------------------------------------------------------------------------
530  '''
531  def __init__(self, Year, Mon,Day, Hr, Minute, Sec, dUT1,lsdat,TimeZone=0,TypeUTIn='UTC'):
532 
533  self.year = Year
534  self.mon = Mon
535  self.day = Day
536  self.hr = Hr
537  self.minute = Minute
538  self.sec = Sec
539  self.typeUTIn = TypeUTIn
540  self.timezone = TimeZone
541  self.localhr = TimeZone + Hr
542  if TypeUTIn == 'UT1' :
543  self.ut1 = hms2ut(self.localhr,Minute,Sec)
544  self.jdut1 = JDay(Year,Mon,Day, self.localhr, Minute, Sec)
545  self.tut1 = (self.jdut1 - 2451545.0 )/ 36525.0
546  self.utc = self.ut1 - dUT1
547  else:
548  self.utc = hms2ut(self.localhr,Minute,Sec)
549  self.ut1 = self.utc + dUT1
550  hrTemp,minTemp,secTemp = ut2hms(self.ut1)
551  self.jdut1 = JDay(Year,Mon,Day, np.int(hrTemp), np.int(minTemp),secTemp)
552  self.tut1 = (self.jdut1 - 2451545.0 )/ 36525.0
553 
554  self.tai = self.utc + lsdat
555  self.tt = self.tai + 32.184
556  hrTemp,minTemp,secTemp = ut2hms(self.tt)
557  self.jdtt = JDay(Year,Mon,Day, hrTemp, minTemp, secTemp)
558  self.ttt = (self.jdtt - 2451545.0 )/ 36525.0
559 
560  #self.me = 357.5277233 + 35999.05034*self.ttt # approx - should do with TTDB
561  self.me = np.fmod((357.5277233 + 35999.05034*self.ttt),360.0)*deg2Rad
562  self.tdb = self.tt + 0.001658 * np.sin(self.me) + 0.00001385*np.sin(2.0*self.me)
563  hrTemp,minTemp,secTemp = ut2hms(self.tdb)
564  self.jdtdb = JDay(Year,Mon,Day, hrTemp, minTemp, secTemp)
565  self.ttdb = (self.jdtdb - 2451545.0)/ 36525.0
566 
567 def precession(TTT):
568  '''
569 # * ----------------------------------------------------------------------------
570 # *
571 # * SUBROUTINE PRECESSION
572 # *
573 # * this function calulates the transformation matrix for precession.
574 # *
575 # * Author : David Vallado 719-573-2600 1 Mar 2001
576 # *
577 # * Inputs Description Range / Units
578 # * TTT - Julian Centuries of TT centuries
579 # *
580 # * Outputs :
581 # * Prec - Precession transformation (eci-mod)
582 # *
583 # * Locals :
584 # * TTT2 - TTT squared
585 # * TTT3 - TTT cubed
586 # * Zeta - PRECESSION ANGLE rad
587 # * z - PRECESSION ANGLE rad
588 # * Theta - PRECESSION ANGLE rad
589 # *
590 # * Coupling :
591 # * none
592 # *
593 # * References :
594 # * Vallado 2007, 228
595 # *
596 # * ----------------------------------------------------------------------------
597  '''
598 # MJM note: Reid Reynolds uses TTDB not TTT.
599 
600  # --------------------- PRECESSION angles ---------------------
601  # MJM These are in arcsec; look like values in Reid Reynolds memo
602  zeta= (( 0.017998*TTT + 0.30188)*TTT + 2306.2181)*TTT
603  theta = (( - 0.041833*TTT - 0.42665)*TTT + 2004.3109)*TTT
604  z = (( 0.018203*TTT + 1.09468)*TTT + 2306.2181)*TTT
605 
606  zeta = zeta * deg2Rad / 3600.0 # convert to radians
607  theta= theta * deg2Rad / 3600.0
608  z = z * deg2Rad / 3600.0
609 
610  coszeta = np.cos(zeta)
611  sinzeta = np.sin(zeta)
612  costheta = np.cos(theta)
613  sintheta = np.sin(theta)
614  cosz = np.cos(z)
615  sinz = np.sin(z)
616 
617  prec=np.zeros((3,3))
618  # ----------------- form matrix J2000 to MOD -----------------
619  prec[0,0] = coszeta * costheta * cosz - sinzeta * sinz
620  prec[0,1] = -sinzeta * costheta * cosz - coszeta * sinz
621  prec[0,2] = -sintheta * cosz
622 
623  prec[1,0] = coszeta * costheta * sinz + sinzeta * cosz
624  prec[1,1] = -sinzeta * costheta * sinz + coszeta * cosz
625  prec[1,2] = -sintheta * sinz
626 
627  prec[2,0] = coszeta * sintheta
628  prec[2,1] = -sinzeta * sintheta
629  prec[2,2] = costheta
630 
631  return prec
632 
633 def nutation(TTT,nut_data):
634  '''
635  # * ----------------------------------------------------------------------------
636  # *
637  # * SUBROUTINE NUTATION
638  # *
639  # * this function calulates the transformation matrix for nutation.
640  # *
641  # * Author : David Vallado 719-573-2600 1 Mar 2001
642  # *
643  # * Inputs Description Range / Units
644  # * TTT - Julian Centuries of TT
645  # *
646  # * Outputs :
647  # * Nut - Nutation transformation (mod-tod)
648  # * DeltaPsi - NUTATION ANGLE rad
649  # * TrueEps - True obliquity of the ecliptic rad
650  # * Omega - rad
651  # * MeanEps - Mean obliquity of the ecliptic rad
652  # *
653  # * Locals :
654  # * TTT2 - TTT squared
655  # * TTT3 - TTT cubed
656  # * MeanEps - Mean obliquity of the ecliptic rad
657  # * l - rad
658  # * ll - rad
659  # * F - rad
660  # * D - rad
661  # * DeltaEps - Change in obliquity rad
662  # *
663  # * Coupling :
664  # * none
665  # *
666  # * References :
667  # * Vallado 2007, 228
668  # *
669  # * ----------------------------------------------------------------------------
670  '''
671  # ---- Determine coefficients for IAU 1980 NUTATION Theory ----
672  TTT2= TTT*TTT
673  TTT3= TTT2*TTT
674  TTT4= TTT2*TTT2
675 
676 # Meaneps is in arcseconds first
677 
678  MeanEps = ((0.001813*TTT - 0.00059)*TTT - 46.8150)*TTT + \
679  84381.448
680 
681  MeanEps = np.fmod(MeanEps/3600.0 , 360.0) # degrees, [0-360)
682  MeanEps = MeanEps * deg2Rad # radians
683 
684  l = 134.96340251 + ( 1717915923.2178*TTT + \
685  31.8792*TTT2 + 0.051635*TTT3 - 0.00024470*TTT4 ) \
686  / 3600.0
687  l1 = 357.52910918 + ( 129596581.0481*TTT - \
688  0.5532*TTT2 - 0.000136*TTT3 - 0.00001149*TTT4 ) \
689  / 3600.0
690  F = 93.27209062 + ( 1739527262.8478*TTT - \
691  12.7512*TTT2 + 0.001037*TTT3 + 0.00000417*TTT4 ) \
692  / 3600.0
693  D = 297.85019547 + ( 1602961601.2090*TTT - \
694  6.3706*TTT2 + 0.006593*TTT3 - 0.00003169*TTT4 ) \
695  / 3600.0
696  Omega= 125.04455501 + ( -6962890.2665*TTT + \
697  7.4722*TTT2 + 0.007702*TTT3 - 0.00005939*TTT4 ) \
698  / 3600.0
699 # ! Above are in degrees and they correspond to:
700 # ! l=mean anomaly of moon
701 # ! l1 = mean anomaly of sun
702 # ! F = mean latitude of moon
703 # ! D = mean elongation of sun
704 # ! omega = ascending node of the moon
705 # ! Because the coefficients (of above in degrees) are all integers,
706 # ! we can safely MOD to get into range of 0-360; but not necessary.
707  l = np.fmod( l , 360.0 ) * deg2Rad
708  l1 = np.fmod( l1 , 360.0 ) * deg2Rad
709  F = np.fmod( F , 360.0 ) * deg2Rad
710  D = np.fmod( D , 360.0 ) * deg2Rad
711  Omega= np.fmod( Omega , 360.0 ) * deg2Rad
712 
713  DeltaPsi= 0.0
714  DeltaEps= 0.0
715 
716  for i in range(len(nut_data)-1,-1,-1):
717  TempVal= nut_data['a1'][i]*l + nut_data['a2'][i]*l1 + nut_data['a3'][i]*F + \
718  nut_data['a4'][i]*D + nut_data['a5'][i]*Omega
719  DeltaPsi= DeltaPsi + (nut_data['A'][i]+ nut_data['B'][i]*TTT) * \
720  np.sin( TempVal )
721  DeltaEps= DeltaEps + (nut_data['C'][i]+ nut_data['D'][i]*TTT) * \
722  np.cos( TempVal )
723 
724 # --------------- Find NUTATION Parameters --------------------
725 # ! MJM begin
726 # ! Now, Reid Reynolds document with the SAME values says they are in
727 # ! microarcseconds and microarcseconds per century (RAR80 terms).
728 # ! other references show these units (for same values!!!) as 0.1mas
729 # ! deltaPsi = deltaPsi* 1.d-4 / 3600. ! convert to as, then deg
730 # ! deltaEps = deltaEps* 1.d-4 / 3600. ! convert to as, then deg
731 # ! ACTUALLY THESE ARE CONVERTED WHEN THEY ARE READ!!!! See INITREDUC
732 # ! MJM end
733  DeltaPsi = np.fmod(DeltaPsi , 360.0 ) * deg2Rad
734  DeltaEps = np.fmod(DeltaEps , 360.0 ) * deg2Rad
735  TrueEps = MeanEps + DeltaEps
736 
737 # Checked order against Reid Reynolds. Multiplication is correct.
738  cospsi = np.cos(DeltaPsi)
739  sinpsi = np.sin(DeltaPsi)
740  coseps = np.cos(MeanEps)
741  sineps = np.sin(MeanEps)
742  costrueeps = np.cos(TrueEps)
743  sintrueeps = np.sin(TrueEps)
744 
745  nut = np.zeros((3,3))
746  nut[0,0] = cospsi
747  nut[0,1] = -coseps * sinpsi
748  nut[0,2] = -sineps * sinpsi
749 
750  nut[1,0] = costrueeps * sinpsi
751  nut[1,1] = costrueeps * coseps * cospsi + sintrueeps * sineps
752  nut[1,2] = costrueeps * sineps * cospsi - sintrueeps * coseps
753 
754  nut[2,0] = sintrueeps * sinpsi
755  nut[2,1] = sintrueeps * coseps * cospsi - sineps * costrueeps
756  nut[2,2] = sintrueeps * sineps * cospsi + costrueeps * coseps
757 
758  return DeltaPsi, TrueEps, MeanEps, Omega, nut
759 
760 def sidereal(jdut1,DeltaPsi,MeanEps,Omega,LOD,terms):
761  # * ----------------------------------------------------------------------------
762  # *
763  # * SUBROUTINE SIDEREAL
764  # *
765  # * this function calulates the transformation matrix that accounts for the
766  # * effects of nutation. Notice that deltaspi should not be moded to a
767  # * positive number because it is multiplied rather than used in a
768  # * trigonometric argument.
769  # *
770  # * Author : David Vallado 719-573-2600 1 Mar 2001
771  # *
772  # * Inputs Description Range / Units
773  # * JDUT1 - Julian Date of UT1 days from 4713 BC
774  # * DeltaPsi - NUTATION ANGLE rad
775  # * MeanEps - Mean obliquity of the ecliptic rad
776  # * Omega - rad
777  # * LOD - Excess length of day sec
778  # * terms - number of terms to include with ast 0, 2
779  # *
780  # * Outputs :
781  # * St - Sidereal Time transformation (tod-pef)
782  # * StDot - Sidereal Time rate transformation (tod-pef)
783  # * Omegaearth - rotation of the earth rad
784  # *
785  # * Locals :
786  # * GST - Mean Greenwich SIDEREAL Time 0 to 2Pi rad
787  # * AST - Apparent GST 0 to 2Pi rad
788  # * Hr - hour hr
789  # * minute - minutes minute
790  # * SEC - seconds SEC
791  # *
792  # * Coupling :
793  # * none
794  # *
795  # * References :
796  # * Vallado 2007, 228
797  # *
798  # * ----------------------------------------------------------------------------
799 
800  OmegaEarth = np.zeros(3)
801 
802  Conv1 = np.pi / (180.0*3600.0)
803 
804  # ------------------------ Find Mean GST ----------------------
805  gst= gstime( jdut1 )
806 
807  # ------------------------ Find Mean AST ----------------------
808  if terms > 0 and jdut1 > 2450449.5:
809  ast = gst + DeltaPsi* np.cos(MeanEps) \
810  + 0.00264*Conv1*np.sin(Omega) \
811  + 0.000063*Conv1*np.sin(2.0*Omega)
812  else:
813  ast = gst + DeltaPsi* np.cos(MeanEps)
814 
815  st = np.zeros((3,3))
816  st[0,0] = np.cos(ast)
817  st[0,1] = np.sin(ast)
818  st[0,2] = 0.
819 
820  st[1,0] = -np.sin(ast)
821  st[1,1] = np.cos(ast)
822  st[1,2] = 0.
823 
824  st[2,0] = 0.
825  st[2,1] = 0.
826  st[2,2] = 1.
827 
828  # ------------ compute sidereal time rate matrix --------------
829  ThetaSa = 7.29211514670698e-05 * (1.0 - LOD/86400.0)
830  OmegaEarth[2] = ThetaSa
831 
832  stdot = np.zeros((3,3))
833  stdot[0,0] = -OmegaEarth[2] * np.sin(ast)
834  stdot[0,1] = OmegaEarth[2] * np.cos(ast)
835  stdot[0,2] = 0.0
836 
837  stdot[1,0] = -OmegaEarth[2] * np.cos(ast)
838  stdot[1,1] = -OmegaEarth[2] * np.sin(ast)
839  stdot[1,2] = 0.0
840 
841  stdot[2,0] = 0.0
842  stdot[2,1] = 0.0
843  stdot[2,2] = 0.0
844 
845  return st,stdot,OmegaEarth
846 
847 def gstime(jd):
848  # * -----------------------------------------------------------------------------
849  # *
850  # * FUNCTION GSTIME
851  # *
852  # * this function finds the Greenwich sidereal time (iau-82).
853  # *
854  # * Author : David Vallado 719-573-2600 1 Mar 2001
855  # *
856  # * Inputs Description Range / Units
857  # * JD - Julian Date days from 4713 BC
858  # *
859  # * OutPuts :
860  # * GSTIME - Greenwich SIDEREAL Time 0 to 2Pi rad
861  # *
862  # * Locals :
863  # * Temp - Temporary variable for reals rad
864  # * TUT1 - Julian Centuries from the
865  # * Jan 1, 2000 12 h epoch (UT1)
866  # *
867  # * Coupling :
868  # *
869  # * References :
870  # * vallado 2007, 193, Eq 3-43
871  # * -----------------------------------------------------------------------------
872 
873 
874  TUT1= ( jd - 2451545.0 ) / 36525.0
875 # !!! MJM Commented out this implmenetation !!!!!
876 # T2= - 6.2D-6*TUT1*TUT1*TUT1
877 # & + 0.093104D0*TUT1*TUT1
878 # & + (876600.0D0*3600.0D0 + 8640184.812866D0)*TUT1
879 # & + 67310.54841D0
880 # T2= DMOD( T2*Deg2Rad/240.0D0,TwoPi ) ! 360/86400 = 1/240, to deg, to rad
881 # !!! END MJM comment out
882  TUT1=0.13090052920192846
883 # Begin MJM add (based on Reid Reynolds), better numerically
884  temp= 280.46061837 + 360.98564736629 * (jd - 2451545.0) + \
885  ( 0.000387933 - (TUT1/38710000.0))*TUT1**2
886  temp=np.fmod(temp*deg2Rad,2*np.pi )
887 # END MJM added
888 
889  # ------------------------ Check quadrants --------------------
890  if temp < 0.0:
891  temp += 2*np.pi
892 
893  return temp
894 
895 def ecef2latlon(xyz):
896  # ! Need to check this set of equations; it seems pretty good,
897  # ! but now there are different equations on Wikipedia. See the papers
898  # ! I downloaded to the ReferenceFrames folder.
899  #
900  # ! MJM copied algorithm from web page
901  # ! http://en.wikipedia.org/wiki/Geodetic_system
902  # ! COnverted to IDL 2009 Sept. 29.
903  # ! Converted to F90 2012 Oct 25
904  # ! 2012/10/31 Converted to Vermeille's algorithm since at least I have a referece. But performance looks
905  # ! exactly the same.
906  # ! 2012/10/31 Changed output order to LON (deg), LAT (DEG), height (m)
907  # !
908  # ! ASSUMES WGS 84 ellipsoid
909  # ! Results are GEODETIC latitude and longitude
910  # ! These will be different from Geocentric latitude
911  a=6378137.0
912  a2=a**2
913  f=1./298.257223563
914  e2= 2*f-f**2
915  e=np.sqrt(e2)
916  e4=e2**2
917 # ! Vermeille's Algorithm (Journal of Geodesy (2002) 76:451-454
918 # ! e=sqrt(e2)
919 # ! a2=a**2
920 
921  p=(xyz[0,:]**2 + xyz[1,:]**2)/a2
922  q=(1-e2)/a2*xyz[2,:]**2
923  r=(p+q-e4)/6
924  s=e4*p*q/4/r**3
925  t=(1. + s + np.sqrt(s*(2. + s)))**(1./3.)
926  u=r*(1. + t + 1./t)
927  v=np.sqrt(u**2+q*e4)
928  w=e2*(u+v-q)/2/v
929  k=np.sqrt(u+v+w**2)-w
930  d=k*np.sqrt(xyz[0,:]**2 + xyz[1,:]**2)/(k+e2)
931 
932  tmp=np.sqrt(d**2 + xyz[2,:]**2)
933 
934  llh = np.zeros((3,len(xyz[0,:])))
935  llh[0,:] = np.arctan2(xyz[1,:],xyz[0,:]) # longitude in rad
936  llh[1,:] = 2*np.arctan2(xyz[2,:],(d + tmp)) # latitude in rad
937  llh[2,:] = tmp*(k+e2-1)/k # height in same units as "a", i.e., meters
938 
939  return llh
940 
941 def wgs84_intercept(rsc,u):
942  # rsc ! spacecraft radius
943  # u ! unit vector pointing toward earth, a lot of them
944  ff=1.0/298.257223563 # WGS-84
945  re=6378137. # WGS-84
946  # ! F is really the diagonal elements of a 3x3 matrix, all rest are 0
947  # ! but we only need the diagonal for the WGS-84 ellipsoid fo rvery simple
948  # ! mathematics, below.
949  F=[1.,1.,1./(1.-ff)**2]
950 
951  # ! Using F as above (but see the comments - really F as the 3x3 matrix)
952  # ! If Rg is a vector from the center to the surface of the earth, then
953  # ! the equation of the WGS-84 ellipsoid is transpose(Rg)F*Rg=Re*Re
954  #
955  # ! The equation of a ray with unit direction (vector) u from the spacecraft
956  # ! (at vector position Rsc) to the earth is the vector equation:
957  # ! Rg = Rsc + s*u ; s is length of ray from the spacecraft to earth;
958  # ! substitue this vector equation into the one above, and use the form
959  # ! of the various vectors and matrices to get the equations below,
960  # ! which yields a simple quadratic equation for s. For us, Rsc is one
961  # ! location, and we have 512 unit direction vectors (elements of u).
962  #
963  # ! precalculate a few things, using the magic rules of sum () and matmul()
964 
965  c=sum(F*rsc**2)-re**2 # scalar, POSITIVE since spacecraft is above earth surface
966 
967  b=2.*np.dot(F*rsc,u) # ns elements; negative since we're looking at earth and this is
968 
969  # ! essentially a dot prod of spacecraft vector (away from center of earth) and view vector
970  # ! (towards earth).
971 
972  a=np.dot(F,np.square(u)) # ns elements, positive, since it is like u^2.
973 
974  det= np.square(b) - 4*np.dot(a,c)
975 
976  if (det < 0).any():
977  print('ERROR IN WGS84_intercept: invalid answer. no intercept')
978  print('CHECK INPUT!')
979  print((det < 0))
980  sys.exit
981 
982 
983  # ! Note that -b is positive. So the closest root, the one with the smallest s,
984  # ! that is the smallest distance from the spacecraft (the one on the spacecraft
985  # ! side of the earth) is the one with the negative sign before the sqrt(). The OTHER one (positive)
986  # ! is when the ray emerges from the earth.
987  # ! Now the distance along the ray, from the spacecraft to the near intercept is:
988  s= (-b -np.sqrt(det))/(2. * a ) # ns elements
989  # ! Once you know s, rout is the the location in ECEF of the intercept with the Earth,
990  # ! traveling a distance s from rsc in the direction(s) u.
991  nr = len(u[:,0])
992  ns = len(u[0,:])
993  rout = np.zeros((nr,ns))
994  #rout= spread(rsc,dim=2,ncopies=ns) + spread(s,dim=1,ncopies=3)*u # 3xns
995  sout=np.tile(s[:,np.newaxis],(1,nr)).T*u
996  rscout = np.tile(rsc[np.newaxis,:],(ns,1)).T
997  rout=np.tile(rsc[np.newaxis,:],(ns,1)).T + sout #*u
998  # Now we need an ECEF->LATLONH conversion
999  out = ecef2latlon(rout)
1000 
1001  # ! The below essentially use terms from the ECEF -> ENU conversion on the surface of an oblate spheroid,
1002  # ! our WGS-84 ellipsoid
1003  # ! The normal to the ellipsoid has the normal snorm
1004  snorm = np.zeros((nr,ns))
1005  snorm[0,:] =np.cos(out[0,:])*np.cos(out[1,:])
1006  snorm[1,:] =np.sin(out[0,:])*np.cos(out[1,:])
1007  snorm[2,:] =np.sin(out[1,:])
1008  # The cos(view zenith) is formed by the pointing vector from the ground to the spacecraft
1009  # dotted into the surface normal; this is (for one location) sum(-u(:,i)*snorm(:,i))
1010  view_zen=np.arccos(np.sum(-u*snorm,axis=0))*rad2Deg # deg from zenith
1011 
1012  uDotE=np.sin(out[0,:])*u[0,:] - np.cos(out[0,:])*u[1,:]
1013  uDotN=np.cos(out[0,:])*np.sin(out[1,:])*u[0,:] + np.sin(out[0,:])*np.sin(out[1,:])*u[1,:] - \
1014  np.cos(out[1,:])*u[2,:]
1015  view_az=np.arctan2(uDotE,uDotN)*rad2Deg # deg from N, clockwise
1016  out=out*rad2Deg
1017 
1018  return out,view_zen,view_az
1019 
1020 def solar_geometry(iyr,imn,idy,xh,xm,xs,xlat,xlong):
1021  # ! 2012-10-31 MJM Made able to handle an array of xlat and xlong, where xlat and xlong
1022  # ! are each vectors
1023  # ! 2002-12-11 MJM Turned into array valued function. Now should only
1024  # ! necessary solar geometry calculations. This is
1025  # ! called many nlines times when given appropriate info for
1026  # ! each line.
1027  # ! xlat: positive NORTH
1028  # ! xlon: positive WEST
1029  # ! 2001-05-16 MJM Solar factors only.
1030 
1031  sol_zen_az = np.zeros((2,len(xlat)))
1032 
1033 
1034  tt=2*np.pi*((xh)/24.+xm/1440.+xs/86400.)
1035 
1036  xlatr = xlat * deg2Rad
1037  xlongr = xlong * deg2Rad
1038 
1039  # at this TIME there is only ONE astronomical solar position
1040  dec,haz = suncor(idy,imn,iyr,tt)
1041  solaz,el = hazel(haz+tt-xlongr,dec,xlatr)
1042 
1043  #!!!C---Note: DEC, SOLAZ,and EL DEC are in radians
1044 
1045  # ! IF (EL .LE. 0.) THEN
1046  # ! print*,IYR,Imn,idy,IH,IM,IS,xlat,xlong
1047  # ! print*,'el=',el
1048  # ! call fatal_error(file_name,module_name,subroutine_name,&
1049  # ! &'ERROR: Sun is below the horizon!!!'//&
1050  # ! &' Check input date, time, latitude and longitude.')
1051  # ! ENDIF
1052 
1053  solzni = np.pi/2.0 - el
1054  sol_zen_az[0,:]=solzni #! radians
1055  sol_zen_az[1,:]=solaz #!radians
1056 
1057  return sol_zen_az
1058 
1059 def suncor(iday,month,iyr,tz):
1060 
1061  jd=julian(iday,month,iyr)
1062  fjd=0.5 + tz/(2.0*np.pi)
1063  ras,dec,gsdt,bzero,pzero,solong = solcor(jd,fjd)
1064  haz=gsdt-ras-tz
1065 
1066  return dec,haz
1067 
1068 def julian(iday,month,iyr):
1069 
1070  md=[0,31,59,90,120,151,181,212,243,273,304,334]
1071 
1072  jyr=iyr-1600
1073  I1=np.int(jyr/400) #Number of complete 400 yr spans; each have 97 leap years
1074  I2=np.int((jyr-400*I1)/100) #Number of 100 year spans, taking away the 400 year spans above
1075  I3=np.int((jyr-400*I1-100*I2)/4) #Number of four year spans, taking away complete 400 and 100 above
1076  # JD(AD1600) normal yr 400 yr 100yr leap yr this century
1077  jday=2305447 + 365*jyr + 97*I1 + 24*I2 +I3 #Counts IYR if IYR a leap year
1078 
1079  jday=jday+md[month-1]+iday
1080 
1081  if month == 2:
1082  leap = 0
1083  if isLeapYear(jyr):
1084  leap = 1
1085 
1086  #Already counted leap day above, so take out if jan/feb
1087  jday = jday - leap
1088  return jday
1089 
1090 def hazel(h,d,xlat):
1091  sne = np.sin(d)*np.sin(xlat)+np.cos(d)*np.cos(xlat)*np.cos(h)
1092  e=np.arcsin(sne)
1093  sna = np.cos(d)*np.sin(h)
1094  csa=(np.sin(xlat)*np.cos(h)*np.cos(d)-np.sin(d)*np.cos(xlat))
1095  a=np.arctan2(sna,csa)+np.pi
1096 
1097  return a,e
1098 
1099 def solcor(jd,fjd):
1100 
1101  jyr=365.25
1102  d=(jd-2415020)+fjd
1103  iyr=np.int(d/jyr)
1104 
1105  g=-.026601523+.01720196977*d-1.95e-15*d*d -2*np.pi*iyr
1106  xlms=4.881627938+.017202791266*d+3.95e-15*d*d-2*np.pi*iyr
1107  obl=.409319747-6.2179e-9*d
1108  ecc=.01675104-1.1444e-9*d
1109  f=d-jyr*iyr
1110  dgsdt=1.739935476+2*np.pi*f/jyr+1.342027e-4*d/jyr
1111  gsdt=dgsdt+2*np.pi*(fjd-0.5)
1112  xlts=xlms+2.*ecc*np.sin(g)+1.25*ecc*ecc*np.sin(2.*g)
1113  sndc=np.sin(xlts)*np.sin(obl)
1114  ddecs=np.arcsin(sndc)
1115  decs=ddecs
1116  csra=np.cos(xlts)/np.cos(ddecs)
1117  ras=np.arccos(csra)
1118  if np.sin(xlts) < 0.:
1119  ras=2*np.pi-ras
1120  omega=1.297906+6.66992e-7*d
1121  thetac=xlts-omega
1122  bzro=np.arcsin(.126199*np.sin(thetac))
1123  p=-np.arctan(np.cos(xlts)*np.tan(obl))-np.arctan(.127216*np.cos(thetac))
1124  xlmm=np.arctan2(.992005*np.sin(thetac),np.cos(thetac))
1125  jdr=jd-2398220
1126  frot=(jdr+fjd)/25.38-np.int((jdr+fjd)/25.38)
1127  solong=xlmm-2*np.pi*frot+np.pi-3.e-4
1128  if solong < 0.:
1129  solong += 2*np.pi
1130  elif solong > 2*np.pi:
1131  solong -= 2*np.pi
1132 
1133  return ras,decs,gsdt,bzro,p,solong
1134 
1135 def frac2DegMinSec(fracd,fracm,fracs):
1136 # ! A subroutine to get degrees, min, sec from
1137 # ! fractional pieces. The output are all real, although a
1138 # ! strong case can be made for degrees and minutes to be
1139 # ! integers.
1140 # ! MJM Unknown date
1141  d_all = fracd+fracm/60.+fracs/3600.
1142  d_hold = np.floor(fracd+fracm/60.+fracs/3600.)
1143  m_all = (d_all-d_hold)*60.
1144  m_hold = np.floor(m_all)
1145 
1146  return [d_hold,m_hold,(m_all-m_hold)*60.]
hico.astreduc.initReduc
def initReduc(filename)
Definition: astreduc.py:188
hico.astreduc.ut2hms
def ut2hms(UTSec)
Definition: astreduc.py:384
hico.astreduc.sidereal
def sidereal(jdut1, DeltaPsi, MeanEps, Omega, LOD, terms)
Definition: astreduc.py:760
hico.astreduc.getEDFData
def getEDFData(filename, year, mon, day)
Definition: astreduc.py:200
hico.astreduc.suncor
def suncor(iday, month, iyr, tz)
Definition: astreduc.py:1059
hico.astreduc.gstime
def gstime(jd)
Definition: astreduc.py:847
hico.astreduc.nutation
def nutation(TTT, nut_data)
Definition: astreduc.py:633
hico.astreduc.solar_geometry
def solar_geometry(iyr, imn, idy, xh, xm, xs, xlat, xlong)
Definition: astreduc.py:1020
hico.astreduc.hms2ut
def hms2ut(Hr, Minute, Sec)
Definition: astreduc.py:348
hico.astreduc.wgs84_intercept
def wgs84_intercept(rsc, u)
Definition: astreduc.py:941
hico.astreduc.parseLSfile
def parseLSfile(filename, skiphdr=0)
Definition: astreduc.py:216
hico.astreduc.frac2DegMinSec
def frac2DegMinSec(fracd, fracm, fracs)
Definition: astreduc.py:1135
hico.astreduc.increment_seconds_to_date
def increment_seconds_to_date(delta)
Definition: astreduc.py:131
hico.astreduc.precession
def precession(TTT)
Definition: astreduc.py:567
hico.astreduc.UT2time.__init__
def __init__(self, Year, Mon, Day, Hr, Minute, Sec, dUT1, lsdat, TimeZone=0, TypeUTIn='UTC')
Definition: astreduc.py:531
hico.astreduc.getLSData
def getLSData(filename, year, mon, day)
Definition: astreduc.py:210
hico.astreduc.julian
def julian(iday, month, iyr)
Definition: astreduc.py:1068
hico.astreduc.parseEDFfile
def parseEDFfile(filename, skiphdr=0)
Definition: astreduc.py:203
hico.astreduc.hazel
def hazel(h, d, xlat)
Definition: astreduc.py:1090
hico.astreduc.polarm
def polarm(xp, yp)
Definition: astreduc.py:254
hico.astreduc.quat_to_rot
def quat_to_rot(q)
Definition: astreduc.py:278
hico.astreduc.solcor
def solcor(jd, fjd)
Definition: astreduc.py:1099
hico.astreduc.ecef2latlon
def ecef2latlon(xyz)
Definition: astreduc.py:895
hico.astreduc.isLeapYear
def isLeapYear(yyyy)
Definition: astreduc.py:117
color_dtdb.range
range
Definition: color_dtdb.py:331
hico.astreduc.JDay
def JDay(Year, Mon, Day, Hr, minute, Sec)
Definition: astreduc.py:425
hico.astreduc.get_astrodate
def get_astrodate(delta_seconds)
Definition: astreduc.py:112


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Level 3 Time Series Plotter

Request a SeaHawk Image

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster