2 Created on Nov 24, 2015
4 @author: rhealy (richard.healy@nasa.gov)
10 from datetime
import datetime
as DT
11 from datetime
import timedelta
as TD
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}
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']
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']
29 as2Rad = as2Deg*deg2Rad
31 if __name__ == '__main__':
34 ocvarroot = os.environ['OCVARROOT']
36 print("OCSSW environement variable OCVARROOT not set.")
37 print("Typically it is set to $OCSSW/run/var.")
39 filein = "{}/hico/{}".format(ocvarroot,'nutation.dat')
41 nut = initReduc(filein)
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])
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
59 self.hour= np.dtype('int32')
60 self.min = np.dtype('int32')
64 if idinc > ls_2012_06_30:
65 #date after leap second, so adjust this to new epoch
69 idinc=idinc-ls_2012_06_30
70 dinc=dinc-ls_2012_06_30
77 jdinc=np.dtype('int32')
78 while idinc >= 86400*mmlen:
80 if isLeapYear(self.yyyy) == 1:
88 dinc=dinc - 86400*mmlen
96 # at this point, we have less a full month. Let's figure out how many whole days
97 nd=np.int(idinc/86400)
99 idinc=idinc - 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
113 d_epoch = DT(year=2009, month=1, day=1)
114 delta_t = TD(seconds=delta_seconds)
115 return d_epoch + delta_t
119 if (yyyy % 100) == 0:
122 if (yyyy % 400) == 0:
135 mlen = np.array([31,28,31,30,31,30,31,31,30,31,30,31])
137 idinc=
int(np.floor(delta))
139 ls_2012_06_30 = 86400*((365*3)+(31+29+31+30+31+30))+1
144 if idinc > ls_2012_06_30:
149 idinc=idinc-ls_2012_06_30
150 dinc=dinc-ls_2012_06_30
153 jdinc=np.dtype(
'int32')
154 while idinc >= 86400*mmlen:
163 dinc=dinc - 86400*mmlen
173 idinc=idinc - 86400*nd
177 idinc=idinc - hour*3600
178 dinc=dinc - hour*3600
186 return yyyy,mm,dd,hour,minute,sec
190 Convrt= 0.0001e0/3600.0e0
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
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 )])
204 return np.genfromtxt(filename,skip_header=skiphdr,
205 delimiter = edf_delimiter,
206 usecols = edf_usecols,
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] ) )])
217 return np.genfromtxt(filename,skip_header=skiphdr,
218 delimiter = ls_delimiter,
219 usecols = ls_usecols,
224 converted to python by R. Healy 11/27/2015
225 * ----------------------------------------------------------------------------
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.
233 * Author : David Vallado 719-573-2600 1 Mar 2001
235 * Inputs Description Range / Units
236 * xp - Polar motion coefficient rad
237 * yp - Polar motion coefficient rad
240 * PM - Polar motion transformation (pef-ecef)
246 * ROT2 - Rotation about the second axis
247 * ROT3 - Rotation about the third axis
252 * ----------------------------------------------------------------------------
262 pm = np.zeros((3,3),dtype=np.float64)
265 pm[0][1] = sinxp * sinyp
266 pm[0][2] = sinxp * cosyp
273 pm[2][1] = cosxp * sinyp
274 pm[2][2] = cosxp * cosyp
280 #*****************************************************************************80
281 # Converted to Python by R. Healy 11/30/2015
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.
290 # This code is distributed under the GNU LGPL license.
302 # James Foley, Andries van Dam, Steven Feiner, John Hughes,
303 # Computer Graphics, Principles and Practice,
305 # Addison Wesley, 1990.
309 # Input, real ( kind = 8 ) Q(4), the quaternion representing the rotation.
311 # Output, real ( kind = 8 ) a[3,3), the rotation matrix.
314 sin_phi = np.sqrt ( np.sum ( np.square(q[1:4]) ) )
318 angle = 2.0 * np.arctan2( sin_phi, cos_phi )
329 ca = np.cos ( angle )
330 sa = np.sin ( angle )
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
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
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 )
350 * Converted to python by R. Healy 11/30/2015
351 * broken into two functions from fortran sub HMS_SEC:
354 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
358 * this subroutine converts Hours, Minutes and Seconds into seconds from the
359 * beginning of the day.
361 * Author : David Vallado 719 - 573 - 2600 1 Mar 2001
363 * Inputs Description Range / Units
365 * minute - Minutes 0 .. 59
366 * Sec - Seconds 0.0D0 .. 59.99D0
367 * Direction - Which set of vars to output FROM TOO
370 * Sec - Seconds 0.0D0 .. 86400.0D0
373 * Temp - Temporary variable
378 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
380 UTSec= Hr*3600.0 + Minute*60.0 + Sec
386 * Converted to python by R. Healy 11/30/2015
387 * broken into two functions from fortran sub HMS_SEC:
390 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
394 * this subroutine converts Hours, Minutes and Seconds into seconds from the
395 * beginning of the day.
397 * Author : David Vallado 719 - 573 - 2600 1 Mar 2001
399 * Inputs Description Range / Units
401 * minute - Minutes 0 .. 59
402 * Sec - Seconds 0.0D0 .. 59.99D0
403 * Direction - Which set of vars to output FROM TOO
406 * Sec - Seconds 0.0D0 .. 86400.0D0
409 * Temp - Temporary variable
414 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
418 Hr = np.floor( Temp )
419 Minute = np.floor( (Temp - Hr)*60.0 )
420 Sec = (Temp - Hr - Minute/60.0 ) * 3600.0
422 return Hr, Minute, Sec
425 def JDay ( Year,Mon,Day,Hr,minute, Sec ):
427 * Converted to python by R. Healy
428 * -----------------------------------------------------------------------------
432 * this subroutine finds the Julian date given the Year, Month, Day, and Time.
434 * Author : David Vallado 719-573-2600 1 Mar 2001
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'
446 * JD - Julian Date days from 4713 BC
449 * B - Var to aid Gregorian dates
455 * Vallado 2007, 189, Alg 14, Ex 3-14
456 * -----------------------------------------------------------------------------
459 jday = 367.0 * Year \
460 - np.int( (7* (Year+np.int ( (Mon+9)/12) ) ) * 0.25 ) \
461 + np.int( 275*Mon / 9 ) \
463 + ( (Sec/60.0 + minute ) / 60.0 + Hr ) / 24.0
468 * Converted to python by R. Healy 11/30/2015
469 * ------------------------------------------------------------------------------
471 * SUBROUTINE CONVTIME
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.
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
483 * Author : David Vallado 719-573-2600 1 Mar 2001
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
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
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
523 * HMS_SEC - Conversion between hr-minute-SEC .and. seconds
524 * jday - Find the Julian date
527 * vallado 2007, 201, alg 16, ex 3-7
529 * ------------------------------------------------------------------------------
531 def __init__(self, Year, Mon,Day, Hr, Minute, Sec, dUT1,lsdat,TimeZone=0,TypeUTIn='UTC'):
542 if TypeUTIn ==
'UT1' :
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
556 hrTemp,minTemp,secTemp =
ut2hms(self.
tt)
557 self.
jdtt =
JDay(Year,Mon,Day, hrTemp, minTemp, secTemp)
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)
569 # * ----------------------------------------------------------------------------
571 # * SUBROUTINE PRECESSION
573 # * this function calulates the transformation matrix for precession.
575 # * Author : David Vallado 719-573-2600 1 Mar 2001
577 # * Inputs Description Range / Units
578 # * TTT - Julian Centuries of TT centuries
581 # * Prec - Precession transformation (eci-mod)
584 # * TTT2 - TTT squared
586 # * Zeta - PRECESSION ANGLE rad
587 # * z - PRECESSION ANGLE rad
588 # * Theta - PRECESSION ANGLE rad
594 # * Vallado 2007, 228
596 # * ----------------------------------------------------------------------------
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
606 zeta = zeta * deg2Rad / 3600.0
607 theta= theta * deg2Rad / 3600.0
608 z = z * deg2Rad / 3600.0
610 coszeta = np.cos(zeta)
611 sinzeta = np.sin(zeta)
612 costheta = np.cos(theta)
613 sintheta = np.sin(theta)
619 prec[0,0] = coszeta * costheta * cosz - sinzeta * sinz
620 prec[0,1] = -sinzeta * costheta * cosz - coszeta * sinz
621 prec[0,2] = -sintheta * cosz
623 prec[1,0] = coszeta * costheta * sinz + sinzeta * cosz
624 prec[1,1] = -sinzeta * costheta * sinz + coszeta * cosz
625 prec[1,2] = -sintheta * sinz
627 prec[2,0] = coszeta * sintheta
628 prec[2,1] = -sinzeta * sintheta
635 # * ----------------------------------------------------------------------------
637 # * SUBROUTINE NUTATION
639 # * this function calulates the transformation matrix for nutation.
641 # * Author : David Vallado 719-573-2600 1 Mar 2001
643 # * Inputs Description Range / Units
644 # * TTT - Julian Centuries of TT
647 # * Nut - Nutation transformation (mod-tod)
648 # * DeltaPsi - NUTATION ANGLE rad
649 # * TrueEps - True obliquity of the ecliptic rad
651 # * MeanEps - Mean obliquity of the ecliptic rad
654 # * TTT2 - TTT squared
656 # * MeanEps - Mean obliquity of the ecliptic rad
661 # * DeltaEps - Change in obliquity rad
667 # * Vallado 2007, 228
669 # * ----------------------------------------------------------------------------
678 MeanEps = ((0.001813*TTT - 0.00059)*TTT - 46.8150)*TTT + \
681 MeanEps = np.fmod(MeanEps/3600.0 , 360.0)
682 MeanEps = MeanEps * deg2Rad
684 l = 134.96340251 + ( 1717915923.2178*TTT + \
685 31.8792*TTT2 + 0.051635*TTT3 - 0.00024470*TTT4 ) \
687 l1 = 357.52910918 + ( 129596581.0481*TTT - \
688 0.5532*TTT2 - 0.000136*TTT3 - 0.00001149*TTT4 ) \
690 F = 93.27209062 + ( 1739527262.8478*TTT - \
691 12.7512*TTT2 + 0.001037*TTT3 + 0.00000417*TTT4 ) \
693 D = 297.85019547 + ( 1602961601.2090*TTT - \
694 6.3706*TTT2 + 0.006593*TTT3 - 0.00003169*TTT4 ) \
696 Omega= 125.04455501 + ( -6962890.2665*TTT + \
697 7.4722*TTT2 + 0.007702*TTT3 - 0.00005939*TTT4 ) \
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
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) * \
721 DeltaEps= DeltaEps + (nut_data[
'C'][i]+ nut_data[
'D'][i]*TTT) * \
733 DeltaPsi = np.fmod(DeltaPsi , 360.0 ) * deg2Rad
734 DeltaEps = np.fmod(DeltaEps , 360.0 ) * deg2Rad
735 TrueEps = MeanEps + DeltaEps
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)
745 nut = np.zeros((3,3))
747 nut[0,1] = -coseps * sinpsi
748 nut[0,2] = -sineps * sinpsi
750 nut[1,0] = costrueeps * sinpsi
751 nut[1,1] = costrueeps * coseps * cospsi + sintrueeps * sineps
752 nut[1,2] = costrueeps * sineps * cospsi - sintrueeps * coseps
754 nut[2,0] = sintrueeps * sinpsi
755 nut[2,1] = sintrueeps * coseps * cospsi - sineps * costrueeps
756 nut[2,2] = sintrueeps * sineps * cospsi + costrueeps * coseps
758 return DeltaPsi, TrueEps, MeanEps, Omega, nut
760 def sidereal(jdut1,DeltaPsi,MeanEps,Omega,LOD,terms):
800 OmegaEarth = np.zeros(3)
802 Conv1 = np.pi / (180.0*3600.0)
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)
813 ast = gst + DeltaPsi* np.cos(MeanEps)
816 st[0,0] = np.cos(ast)
817 st[0,1] = np.sin(ast)
820 st[1,0] = -np.sin(ast)
821 st[1,1] = np.cos(ast)
829 ThetaSa = 7.29211514670698e-05 * (1.0 - LOD/86400.0)
830 OmegaEarth[2] = ThetaSa
832 stdot = np.zeros((3,3))
833 stdot[0,0] = -OmegaEarth[2] * np.sin(ast)
834 stdot[0,1] = OmegaEarth[2] * np.cos(ast)
837 stdot[1,0] = -OmegaEarth[2] * np.cos(ast)
838 stdot[1,1] = -OmegaEarth[2] * np.sin(ast)
845 return st,stdot,OmegaEarth
874 TUT1= ( jd - 2451545.0 ) / 36525.0
882 TUT1=0.13090052920192846
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 )
921 p=(xyz[0,:]**2 + xyz[1,:]**2)/a2
922 q=(1-e2)/a2*xyz[2,:]**2
925 t=(1. + s + np.sqrt(s*(2. + s)))**(1./3.)
929 k=np.sqrt(u+v+w**2)-w
930 d=k*np.sqrt(xyz[0,:]**2 + xyz[1,:]**2)/(k+e2)
932 tmp=np.sqrt(d**2 + xyz[2,:]**2)
934 llh = np.zeros((3,len(xyz[0,:])))
935 llh[0,:] = np.arctan2(xyz[1,:],xyz[0,:])
936 llh[1,:] = 2*np.arctan2(xyz[2,:],(d + tmp))
937 llh[2,:] = tmp*(k+e2-1)/k
949 F=[1.,1.,1./(1.-ff)**2]
965 c=sum(F*rsc**2)-re**2
972 a=np.dot(F,np.square(u))
974 det= np.square(b) - 4*np.dot(a,c)
977 print(
'ERROR IN WGS84_intercept: invalid answer. no intercept')
978 print(
'CHECK INPUT!')
988 s= (-b -np.sqrt(det))/(2. * a )
993 rout = np.zeros((nr,ns))
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
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,:])
1010 view_zen=np.arccos(np.sum(-u*snorm,axis=0))*rad2Deg
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
1018 return out,view_zen,view_az
1031 sol_zen_az = np.zeros((2,len(xlat)))
1034 tt=2*np.pi*((xh)/24.+xm/1440.+xs/86400.)
1036 xlatr = xlat * deg2Rad
1037 xlongr = xlong * deg2Rad
1040 dec,haz =
suncor(idy,imn,iyr,tt)
1041 solaz,el =
hazel(haz+tt-xlongr,dec,xlatr)
1053 solzni = np.pi/2.0 - el
1054 sol_zen_az[0,:]=solzni
1055 sol_zen_az[1,:]=solaz
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)
1070 md=[0,31,59,90,120,151,181,212,243,273,304,334]
1074 I2=np.int((jyr-400*I1)/100)
1075 I3=np.int((jyr-400*I1-100*I2)/4)
1077 jday=2305447 + 365*jyr + 97*I1 + 24*I2 +I3
1079 jday=jday+md[month-1]+iday
1091 sne = np.sin(d)*np.sin(xlat)+np.cos(d)*np.cos(xlat)*np.cos(h)
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
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
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)
1116 csra=np.cos(xlts)/np.cos(ddecs)
1118 if np.sin(xlts) < 0.:
1120 omega=1.297906+6.66992e-7*d
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))
1126 frot=(jdr+fjd)/25.38-np.int((jdr+fjd)/25.38)
1127 solong=xlmm-2*np.pi*frot+np.pi-3.e-4
1130 elif solong > 2*np.pi:
1133 return ras,decs,gsdt,bzro,p,solong
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)
1146 return [d_hold,m_hold,(m_all-m_hold)*60.]