diff --git a/schainpy/model/graphics/jroplot_base.py b/schainpy/model/graphics/jroplot_base.py index 6562461..42e3ad8 100644 --- a/schainpy/model/graphics/jroplot_base.py +++ b/schainpy/model/graphics/jroplot_base.py @@ -258,6 +258,13 @@ class Plot(Operation): self.tmin = kwargs.get('tmin', None) self.t_units = kwargs.get('t_units', "h_m") self.selectedHeightsList = kwargs.get('selectedHeightsList', []) + self.extFile = kwargs.get('filename', None) + self.bFieldList = kwargs.get('bField', []) + self.celestialList = kwargs.get('celestial', []) + + if isinstance(self.bFieldList, int): + self.bFieldList = [self.bFieldList] + if isinstance(self.selectedHeightsList, int): self.selectedHeightsList = [self.selectedHeightsList] diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index a0b42a4..d6952d1 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -13,6 +13,10 @@ from schainpy.model.graphics.jroplot_base import Plot, plt, log from itertools import combinations from matplotlib.ticker import LinearLocator +from schainpy.model.utils.BField import BField +from scipy.interpolate import splrep +from scipy.interpolate import splev + from matplotlib import __version__ as plt_version if plt_version >='3.3.4': @@ -67,6 +71,11 @@ class SpectraPlot(Plot): norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter noise = 10*numpy.log10(dataOut.getNoise()/norm) +<<<<<<< HEAD +======= + + +>>>>>>> 37cccf17c7b80521b59b978cb30e4ab2e6f37fce z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights)) for ch in range(dataOut.nChannels): if hasattr(dataOut.normFactor,'ndim'): @@ -1466,4 +1475,251 @@ class GeneralProfilePlot(Plot): #self.xmax = max(self.z) ax.plt_r = ax.plot(self.z[i], self.y)[0] else: - ax.plt_r.set_data(self.z[i], self.y) \ No newline at end of file + ax.plt_r.set_data(self.z[i], self.y) + + +########################################################################################################## +########################################## AMISR_V4 ###################################################### + +class RTIMapPlot(Plot): + ''' + Plot for RTI data + + Example: + + controllerObj = Project() + controllerObj.setup(id = '11', name='eej_proc', description=desc) + ##....................................................................................... + ##....................................................................................... + readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader', path=inPath, startDate='2023/05/24',endDate='2023/05/24', + startTime='12:00:00',endTime='12:45:59',walk=1,timezone='lt',margin_days=1,code = code,nCode = nCode, + nBaud = nBaud,nOsamp = nosamp,nChannels=nChannels,nFFT=NFFT, + syncronization=False,shiftChannels=0) + + volts_proc = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId()) + + opObj01 = volts_proc.addOperation(name='Decoder', optype='other') + opObj01.addParameter(name='code', value=code, format='floatlist') + opObj01.addParameter(name='nCode', value=1, format='int') + opObj01.addParameter(name='nBaud', value=nBaud, format='int') + opObj01.addParameter(name='osamp', value=nosamp, format='int') + + opObj12 = volts_proc.addOperation(name='selectHeights', optype='self') + opObj12.addParameter(name='minHei', value='90', format='float') + opObj12.addParameter(name='maxHei', value='150', format='float') + + proc_spc = controllerObj.addProcUnit(datatype='SpectraProc', inputId=volts_proc.getId()) + proc_spc.addParameter(name='nFFTPoints', value='8', format='int') + + opObj11 = proc_spc.addOperation(name='IncohInt', optype='other') + opObj11.addParameter(name='n', value='1', format='int') + + beamMapFile = "/home/japaza/Documents/AMISR_sky_mapper/UMET_beamcodes.csv" + + opObj12 = proc_spc.addOperation(name='RTIMapPlot', optype='external') + opObj12.addParameter(name='selectedHeightsList', value='95, 100, 105, 110 ', format='int') + opObj12.addParameter(name='bField', value='100', format='int') + opObj12.addParameter(name='filename', value=beamMapFile, format='str') + + ''' + + CODE = 'rti_skymap' + + plot_type = 'scatter' + titles = None + colormap = 'jet' + channelList = [] + elevationList = [] + azimuthList = [] + last_noise = None + flag_setIndex = False + heights = [] + dcosx = [] + dcosy = [] + fullDcosy = None + fullDcosy = None + hindex = [] + mapFile = False + ##### BField #### + flagBField = False + dcosxB = [] + dcosyB = [] + Bmarker = ['+','*','D','x','s','>','o','^'] + + + + def setup(self): + + self.xaxis = 'Range (Km)' + if len(self.selectedHeightsList) > 0: + self.nplots = len(self.selectedHeightsList) + else: + self.nplots = 4 + self.ncols = int(numpy.ceil(self.nplots/2)) + self.nrows = int(numpy.ceil(self.nplots/self.ncols)) + self.ylabel = 'dcosy' + self.xlabel = 'dcosx' + self.colorbar = True + self.width = 6 + 4.1*self.nrows + self.height = 3 + 3.5*self.ncols + + + if self.extFile!=None: + try: + pointings = numpy.genfromtxt(self.extFile, delimiter=',') + full_azi = pointings[:,1] + full_elev = pointings[:,2] + self.fullDcosx = numpy.cos(numpy.radians(full_elev))*numpy.sin(numpy.radians(full_azi)) + self.fullDcosy = numpy.cos(numpy.radians(full_elev))*numpy.cos(numpy.radians(full_azi)) + mapFile = True + except Exception as e: + self.extFile = None + print(e) + + + + + def update_list(self,dataOut): + if len(self.channelList) == 0: + self.channelList = dataOut.channelList + if len(self.elevationList) == 0: + self.elevationList = dataOut.elevationList + if len(self.azimuthList) == 0: + self.azimuthList = dataOut.azimuthList + a = numpy.radians(numpy.asarray(self.azimuthList)) + e = numpy.radians(numpy.asarray(self.elevationList)) + self.heights = dataOut.heightList + self.dcosx = numpy.cos(e)*numpy.sin(a) + self.dcosy = numpy.cos(e)*numpy.cos(a) + + if len(self.bFieldList)>0: + datetObj = datetime.datetime.fromtimestamp(dataOut.utctime) + doy = datetObj.timetuple().tm_yday + year = datetObj.year + # self.dcosxB, self.dcosyB + ObjB = BField(year=year,doy=doy,site=2,heights=self.bFieldList) + [dcos, alpha, nlon, nlat] = ObjB.getBField() + + alpha_location = numpy.zeros((nlon,2,len(self.bFieldList))) + for ih in range(len(self.bFieldList)): + alpha_location[:,0,ih] = dcos[:,0,ih,0] + for ilon in numpy.arange(nlon): + myx = (alpha[ilon,:,ih])[::-1] + myy = (dcos[ilon,:,ih,0])[::-1] + tck = splrep(myx,myy,s=0) + mydcosx = splev(ObjB.alpha_i,tck,der=0) + + myx = (alpha[ilon,:,ih])[::-1] + myy = (dcos[ilon,:,ih,1])[::-1] + tck = splrep(myx,myy,s=0) + mydcosy = splev(ObjB.alpha_i,tck,der=0) + alpha_location[ilon,:,ih] = numpy.array([mydcosx, mydcosy]) + self.dcosxB.append(alpha_location[:,0,ih]) + self.dcosyB.append(alpha_location[:,1,ih]) + self.flagBField = True + + if len(self.celestialList)>0: + #getBField(self.bFieldList, date) + #pass = kwargs.get('celestial', []) + pass + + + + def update(self, dataOut): + + if len(self.channelList) == 0: + self.update_list(dataOut) + + if not self.flag_setIndex: + if len(self.selectedHeightsList)>0: + for sel_height in self.selectedHeightsList: + index_list = numpy.where(self.heights >= sel_height) + index_list = index_list[0] + self.hindex.append(index_list[0]) + # else: + # k = len(self.heights) + # self.hindex.append(int(k/2)) + self.flag_setIndex = True + + data = {} + meta = {} + + data['rti_skymap'] = dataOut.getPower() + norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter + noise = 10*numpy.log10(dataOut.getNoise()/norm) + data['noise'] = noise + + return data, meta + + def plot(self): + + ###### + self.x = self.dcosx + self.y = self.dcosy + self.z = self.data[-1]['rti_skymap'] + self.z = numpy.array(self.z, dtype=float) + + #print("inde x1 ", self.height_index) + if len(self.hindex) > 0: + index = self.hindex + else: + index = numpy.arange(0, len(self.heights), int((len(self.heights))/4.2)) + + #print(index) + self.titles = ['Height {:.2f} km '.format(self.heights[i])+" " for i in index] + for n, ax in enumerate(self.axes): + + if ax.firsttime: + + + self.xmax = self.xmax if self.xmax else numpy.nanmax(self.x) + self.xmin = self.xmin if self.xmin else numpy.nanmin(self.x) + + self.ymax = self.ymax if self.ymax else numpy.nanmax(self.y) + self.ymin = self.ymin if self.ymin else numpy.nanmin(self.y) + + self.zmax = self.zmax if self.zmax else numpy.nanmax(self.z) + self.zmin = self.zmin if self.zmin else numpy.nanmin(self.z) + + + if self.extFile!=None: + ax.scatter(self.fullDcosx, self.fullDcosy, marker="+", s=20) + #print(self.fullDcosx) + pass + + + ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin, + s=60, marker="s", vmax = self.zmax) + + + ax.minorticks_on() + ax.grid(which='major', axis='both') + ax.grid(which='minor', axis='x') + + if self.flagBField : + + for ih in range(len(self.bFieldList)): + label = str(self.bFieldList[ih]) + ' km' + ax.plot(self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8], + label=label, linestyle='--', ms=4.0,lw=0.5) + handles, labels = ax.get_legend_handles_labels() + a = -0.05 + b = 1.15 - 1.19*(self.nrows) + self.axes[0].legend(handles,labels, bbox_to_anchor=(a,b), prop={'size': (5.8+ 1.1*self.nplots)}, title='B Field ⊥') + + else: + + ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin, + s=80, marker="s", vmax = self.zmax) + + if self.flagBField : + for ih in range(len(self.bFieldList)): + ax.plot (self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8], + linestyle='--', ms=4.0,lw=0.5) + + # handles, labels = ax.get_legend_handles_labels() + # a = -0.05 + # b = 1.15 - 1.19*(self.nrows) + # self.axes[0].legend(handles,labels, bbox_to_anchor=(a,b), prop={'size': (5.8+ 1.1*self.nplots)}, title='B Field ⊥') + + diff --git a/schainpy/model/graphics/jroplot_voltage.py b/schainpy/model/graphics/jroplot_voltage.py index b6382cd..7c9f82a 100644 --- a/schainpy/model/graphics/jroplot_voltage.py +++ b/schainpy/model/graphics/jroplot_voltage.py @@ -305,6 +305,15 @@ class PulsepairSignalPlot(ScopePlot): class Spectra2DPlot(Plot): ''' Plot for 2D Spectra data + Necessary data as Block + you could use profiles2Block Operation + + Example: + # opObj11 = volts_proc.addOperation(name='profiles2Block', optype='other') + # # opObj11.addParameter(name='n', value=10, format='int') + # opObj11.addParameter(name='timeInterval', value='2', format='int') + + # opObj12 = volts_proc.addOperation(name='Spectra2DPlot', optype='external') ''' CODE = 'spc' diff --git a/schainpy/model/io/jroIO_kamisr.py b/schainpy/model/io/jroIO_kamisr.py index 8d70d40..c45436c 100644 --- a/schainpy/model/io/jroIO_kamisr.py +++ b/schainpy/model/io/jroIO_kamisr.py @@ -61,8 +61,9 @@ class AMISRReader(ProcessingUnit): self.elevationList = [] self.dataShape = None self.flag_old_beams = False - - + + self.flagAsync = False #Use when the experiment has no syncronization + self.shiftChannels = 0 self.profileIndex = 0 @@ -109,6 +110,8 @@ class AMISRReader(ProcessingUnit): ignEndDate=None, ignStartTime=None, ignEndTime=None, + syncronization=True, + shiftChannels=0 ): @@ -123,7 +126,8 @@ class AMISRReader(ProcessingUnit): self.nOsamp = int(nOsamp) self.margin_days = margin_days self.__sampleRate = None - + self.flagAsync = not syncronization + self.shiftChannels = shiftChannels self.nFFT = nFFT self.nChannels = nChannels if ignStartTime!=None and ignEndTime!=None: @@ -178,7 +182,14 @@ class AMISRReader(ProcessingUnit): a = [line for line in linesExp if "nbeamcodes" in line] self.nChannels = int(a[0][11:]) + if not self.flagAsync: #for experiments with no syncronization + self.shiftChannels = 0 + + + self.beamCodeByPulse = fp.get(header+'/BeamCode') # LIST OF BEAMS PER PROFILE, TO BE USED ON REARRANGE + + if (self.startDate > datetime.date(2021, 7, 15)) or self.flag_old_beams: #Se cambió la forma de extracción de Apuntes el 17 o forzar con flag de reorganización self.beamcodeFile = fp['Setup/Beamcodefile'][()].decode() self.trueBeams = self.beamcodeFile.split("\n") @@ -189,10 +200,17 @@ class AMISRReader(ProcessingUnit): beams = [self.trueBeams[b] for b in beams_idx] self.beamCode = [int(x, 16) for x in beams] + if(self.flagAsync and self.shiftChannels == 0): + initBeam = self.beamCodeByPulse[0, 0] + self.shiftChannels = numpy.argwhere(self.beamCode ==initBeam)[0,0] + else: _beamCode= fp.get('Raw11/Data/Beamcodes') #se usa la manera previa al cambio de apuntes self.beamCode = _beamCode[0,:] + + + if self.beamCodeMap == None: self.beamCodeMap = fp['Setup/BeamcodeMap'] for beam in self.beamCode: @@ -219,7 +237,8 @@ class AMISRReader(ProcessingUnit): self.nblocks = self.pulseCount.shape[0] #nblocks self.profPerBlockRAW = self.pulseCount.shape[1] #profiles per block in raw data self.nprofiles = self.pulseCount.shape[1] #nprofile - self.nsa = self.nsamplesPulse[0,0] #ngates + #self.nsa = self.nsamplesPulse[0,0] #ngates + self.nsa = len(self.rangeFromFile[0]) self.nchannels = len(self.beamCode) self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds #print("IPPS secs: ",self.ippSeconds) @@ -521,8 +540,9 @@ class AMISRReader(ProcessingUnit): #profPerCH = int(self.profPerBlockRAW / self.nChannels) for thisChannel in range(nchan): + ich = thisChannel - idx_ch = [self.nFFT*(thisChannel + nchan*k) for k in range(profPerCH)] + idx_ch = [self.nFFT*(ich + nchan*k) for k in range(profPerCH)] #print(idx_ch) if self.nFFT > 1: aux = [numpy.arange(i, i+self.nFFT) for i in idx_ch] @@ -532,14 +552,15 @@ class AMISRReader(ProcessingUnit): else: idx_ch = numpy.array(idx_ch, dtype=int) - #print(thisChannel,profPerCH,idx_ch) - #print(numpy.where(channels==self.beamCode[thisChannel])[0]) - #new_block[:,thisChannel,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[thisChannel])[0],:] - new_block[:,thisChannel,:,:] = self.dataset[:,idx_ch,:] + #print(ich,profPerCH,idx_ch) + #print(numpy.where(channels==self.beamCode[ich])[0]) + #new_block[:,ich,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[ich])[0],:] + new_block[:,ich,:,:] = self.dataset[:,idx_ch,:] new_block = numpy.transpose(new_block, (1,0,2,3)) new_block = numpy.reshape(new_block, (nchan,-1, nsamples)) - + if self.flagAsync: + new_block = numpy.roll(new_block, self.shiftChannels, axis=0) return new_block def updateIndexes(self): @@ -575,10 +596,13 @@ class AMISRReader(ProcessingUnit): # self.dataOut.channelIndexList = None - self.dataOut.azimuthList = numpy.array(self.azimuthList) - self.dataOut.elevationList = numpy.array(self.elevationList) - self.dataOut.codeList = numpy.array(self.beamCode) - + #self.dataOut.azimuthList = numpy.roll( numpy.array(self.azimuthList) ,self.shiftChannels) + #self.dataOut.elevationList = numpy.roll(numpy.array(self.elevationList) ,self.shiftChannels) + #self.dataOut.codeList = numpy.roll(numpy.array(self.beamCode), self.shiftChannels) + + self.dataOut.azimuthList = self.azimuthList + self.dataOut.elevationList = self.elevationList + self.dataOut.codeList = self.beamCode @@ -637,7 +661,11 @@ class AMISRReader(ProcessingUnit): self.dataOut.radarControllerHeaderObj.elevationList = self.elevationList self.dataOut.radarControllerHeaderObj.dtype = "Voltage" self.dataOut.ippSeconds = self.ippSeconds +<<<<<<< HEAD self.dataOut.ippFactor = self.nchannels*self.nFFT +======= + self.dataOut.ippFactor = self.nFFT +>>>>>>> 37cccf17c7b80521b59b978cb30e4ab2e6f37fce pass def readNextFile(self,online=False): diff --git a/schainpy/model/utils/Astro_Coords.py b/schainpy/model/utils/Astro_Coords.py new file mode 100644 index 0000000..5f5554d --- /dev/null +++ b/schainpy/model/utils/Astro_Coords.py @@ -0,0 +1,1420 @@ +""" +The module ASTRO_COORDS.py gathers classes and functions for coordinates transformation. Additiona- +lly a class EquatorialCorrections and celestial bodies are defined. The first of these is to correct +any error in the location of the body and the second to know the location of certain celestial bo- +dies in the sky. + +MODULES CALLED: +OS, NUMPY, NUMERIC, SCIPY, TIME_CONVERSIONS + +MODIFICATION HISTORY: +Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ Sep 20, 2009. +""" + +import numpy +#import Numeric +import scipy.interpolate +import os +import sys +from schainpy.model.utils import TimeTools +from schainpy.model.utils import Misc_Routines + +class EquatorialCorrections(): + def __init__(self): + """ + EquatorialCorrections class creates an object to call methods to correct the loca- + tion of the celestial bodies. + + Modification History + -------------------- + Converted to Object-oriented Programming by Freddy Galindo, ROJ, 27 September 2009. + """ + + pass + + def co_nutate(self,jd,ra,dec): + """ + co_nutate calculates changes in RA and Dec due to nutation of the Earth's rotation + Additionally it returns the obliquity of the ecliptic (eps), nutation in the longi- + tude of the ecliptic (d_psi) and nutation in the pbliquity of the ecliptic (d_eps). + + Parameters + ---------- + jd = Julian Date (Scalar or array). + RA = A scalar o array giving the Right Ascention of interest. + Dec = A scalar o array giving the Right Ascention of interest. + + Return + ------ + d_ra = Correction to ra due to nutation. + d_dec = Correction to dec due to nutation. + + Examples + -------- + >> Julian = 2462088.7 + >> Ra = 41.547213 + >> Dec = 49.348483 + >> [d_ra,d_dec,eps,d_psi,d_eps] = co_nutate(julian,Ra,Dec) + >> print d_ra, d_dec, eps, d_psi, d_eps + [ 15.84276651] [ 6.21641029] [ 0.4090404] [ 14.85990198] [ 2.70408658] + + Modification history + -------------------- + Written by Chris O'Dell, 2002. + Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009. + """ + + jd = numpy.atleast_1d(jd) + ra = numpy.atleast_1d(ra) + dec = numpy.atleast_1d(dec) + + # Useful transformation constants + d2as = numpy.pi/(180.*3600.) + + # Julian centuries from J2000 of jd + T = (jd - 2451545.0)/36525.0 + + # Must calculate obliquity of ecliptic + [d_psi, d_eps] = self.nutate(jd) + d_psi = numpy.atleast_1d(d_psi) + d_eps = numpy.atleast_1d(d_eps) + + eps0 = (23.4392911*3600.) - (46.8150*T) - (0.00059*T**2) + (0.001813*T**3) + # True obliquity of the ecliptic in radians + eps = (eps0 + d_eps)/3600.*Misc_Routines.CoFactors.d2r + + # Useful numbers + ce = numpy.cos(eps) + se = numpy.sin(eps) + + # Convert Ra-Dec to equatorial rectangular coordinates + x = numpy.cos(ra*Misc_Routines.CoFactors.d2r)*numpy.cos(dec*Misc_Routines.CoFactors.d2r) + y = numpy.sin(ra*Misc_Routines.CoFactors.d2r)*numpy.cos(dec*Misc_Routines.CoFactors.d2r) + z = numpy.sin(dec*Misc_Routines.CoFactors.d2r) + + # Apply corrections to each rectangular coordinate + x2 = x - (y*ce + z*se)*d_psi*Misc_Routines.CoFactors.s2r + y2 = y + (x*ce*d_psi - z*d_eps)*Misc_Routines.CoFactors.s2r + z2 = z + (x*se*d_psi + y*d_eps)*Misc_Routines.CoFactors.s2r + + # Convert bask to equatorial spherical coordinates + r = numpy.sqrt(x2**2. + y2**2. + z2**2.) + xyproj =numpy.sqrt(x2**2. + y2**2.) + + ra2 = x2*0.0 + dec2 = x2*0.0 + + xyproj = numpy.atleast_1d(xyproj) + z = numpy.atleast_1d(z) + r = numpy.atleast_1d(r) + x2 = numpy.atleast_1d(x2) + y2 = numpy.atleast_1d(y2) + z2 = numpy.atleast_1d(z2) + ra2 = numpy.atleast_1d(ra2) + dec2 = numpy.atleast_1d(dec2) + + w1 = numpy.where((xyproj==0) & (z!=0)) + w2 = numpy.where(xyproj!=0) + + # Calculate Ra and Dec in radians (later convert to degrees) + if w1[0].size>0: + # Places where xyproj=0 (point at NCP or SCP) + dec2[w1] = numpy.arcsin(z2[w1]/r[w1]) + ra2[w1] = 0 + + if w2[0].size>0: + # Places other than NCP or SCP + ra2[w2] = numpy.arctan2(y2[w2],x2[w2]) + dec2[w2] = numpy.arcsin(z2[w2]/r[w2]) + + # Converting to degree + ra2 = ra2/Misc_Routines.CoFactors.d2r + dec2 = dec2/Misc_Routines.CoFactors.d2r + + w = numpy.where(ra2<0.) + if w[0].size>0: + ra2[w] = ra2[w] + 360. + + # Return changes in Ra and Dec in arcseconds + d_ra = (ra2 -ra)*3600. + d_dec = (dec2 - dec)*3600. + + return d_ra, d_dec, eps, d_psi, d_eps + + def nutate(self,jd): + """ + nutate returns the nutation in longitude and obliquity for a given Julian date. + + Parameters + ---------- + jd = Julian ephemeris date, scalar or vector. + + Return + ------ + nut_long = The nutation in longitude. + nut_obliq = The nutation in latitude. + + Example + ------- + >> julian = 2446895.5 + >> [nut_long,nut_obliq] = nutate(julian) + >> print nut_long, nut_obliq + -3.78793107711 9.44252069864 + + >> julians = 2415020.5 + numpy.arange(50) + >> [nut_long,nut_obliq] = nutate(julians) + + Modification History + -------------------- + Written by W.Landsman (Goddard/HSTX), June 1996. + Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009. + """ + + jd = numpy.atleast_1d(jd) + + # Form time in Julian centuries from 1900 + t = (jd - 2451545.0)/36525.0 + + # Mean elongation of the moon + coeff1 = numpy.array([1/189474.0,-0.0019142,445267.111480,297.85036]) + d = numpy.poly1d(coeff1) + d = d(t)*Misc_Routines.CoFactors.d2r + d = self.cirrange(d,rad=1) + + # Sun's mean elongation + coeff2 = numpy.array([-1./3e5,-0.0001603,35999.050340,357.52772]) + m = numpy.poly1d(coeff2) + m = m(t)*Misc_Routines.CoFactors.d2r + m = self.cirrange(m,rad=1) + + # Moon's mean elongation + coeff3 = numpy.array([1.0/5.625e4,0.0086972,477198.867398,134.96298]) + mprime = numpy.poly1d(coeff3) + mprime = mprime(t)*Misc_Routines.CoFactors.d2r + mprime = self.cirrange(mprime,rad=1) + + # Moon's argument of latitude + coeff4 = numpy.array([-1.0/3.27270e5,-0.0036825,483202.017538,93.27191]) + f = numpy.poly1d(coeff4) + f = f(t)*Misc_Routines.CoFactors.d2r + f = self.cirrange(f,rad=1) + + # Longitude fo the ascending node of the Moon's mean orbit on the ecliptic, measu- + # red from the mean equinox of the date. + coeff5 = numpy.array([1.0/4.5e5,0.0020708,-1934.136261,125.04452]) + omega = numpy.poly1d(coeff5) + omega = omega(t)*Misc_Routines.CoFactors.d2r + omega = self.cirrange(omega,rad=1) + + d_lng = numpy.array([0,-2,0,0,0,0,-2,0,0,-2,-2,-2,0,2,0,2,0,0,-2,0,2,0,0,-2,0,-2,0,0,\ + 2,-2,0,-2,0,0,2,2,0,-2,0,2,2,-2,-2,2,2,0,-2,-2,0,-2,-2,0,-1,-2,1,0,0,-1,0,\ + 0,2,0,2]) + + m_lng = numpy.array([0,0,0,0,1,0,1,0,0,-1]) + m_lng = numpy.append(m_lng,numpy.zeros(17)) + m_lng = numpy.append(m_lng,numpy.array([2,0,2,1,0,-1,0,0,0,1,1,-1,0,0,0,0,0,0,-1,-1,0,0,\ + 0,1,0,0,1,0,0,0,-1,1,-1,-1,0,-1])) + + mp_lng = numpy.array([0,0,0,0,0,1,0,0,1,0,1,0,-1,0,1,-1,-1,1,2,-2,0,2,2,1,0,0, -1, 0,\ + -1,0,0,1,0,2,-1,1,0,1,0,0,1,2,1,-2,0,1,0,0,2,2,0,1,1,0,0,1,-2,1,1,1,-1,3,0]) + + f_lng = numpy.array([0,2,2,0,0,0,2,2,2,2,0,2,2,0,0,2,0,2,0,2,2,2,0,2,2,2,2,0,0,2,0,0,\ + 0,-2,2,2,2,0,2,2,0,2,2,0,0,0,2,0,2,0,2,-2,0,0,0,2,2,0,0,2,2,2,2]) + + om_lng = numpy.array([1,2,2,2,0,0,2,1,2,2,0,1,2,0,1,2,1,1,0,1,2,2,0,2,0,0,1,0,1,2,1, \ + 1,1,0,1,2,2,0,2,1,0,2,1,1,1,0,1,1,1,1,1,0,0,0,0,0,2,0,0,2,2,2,2]) + + sin_lng = numpy.array([-171996,-13187,-2274,2062,1426,712,-517,-386,-301, 217, -158, \ + 129,123,63,63,-59,-58,-51,48,46,-38,-31,29,29,26,-22,21,17,16,-16,-15,-13,\ + -12,11,-10,-8,7,-7,-7,-7,6,6,6,-6,-6,5,-5,-5,-5,4,4,4,-4,-4,-4,3,-3,-3,-3,\ + -3,-3,-3,-3]) + + sdelt = numpy.array([-174.2,-1.6,-0.2,0.2,-3.4,0.1,1.2,-0.4,0,-0.5,0, 0.1, 0, 0, 0.1,\ + 0,-0.1]) + sdelt = numpy.append(sdelt,numpy.zeros(10)) + sdelt = numpy.append(sdelt,numpy.array([-0.1, 0, 0.1])) + sdelt = numpy.append(sdelt,numpy.zeros(33)) + + cos_lng = numpy.array([92025,5736,977,-895,54,-7,224,200,129,-95,0,-70,-53,0,-33,26, \ + 32,27,0,-24,16,13,0,-12,0,0,-10,0,-8,7,9,7,6,0,5,3,-3,0,3,3,0,-3,-3,3,3,0,\ + 3,3,3]) + cos_lng = numpy.append(cos_lng,numpy.zeros(14)) + + cdelt = numpy.array([8.9,-3.1,-0.5,0.5,-0.1,0.0,-0.6,0.0,-0.1,0.3]) + cdelt = numpy.append(cdelt,numpy.zeros(53)) + + # Sum the periodic terms. + n = numpy.size(jd) + nut_long = numpy.zeros(n) + nut_obliq = numpy.zeros(n) + + d_lng = d_lng.reshape(numpy.size(d_lng),1) + d = d.reshape(numpy.size(d),1) + matrix_d_lng = numpy.dot(d_lng,d.transpose()) + + m_lng = m_lng.reshape(numpy.size(m_lng),1) + m = m.reshape(numpy.size(m),1) + matrix_m_lng = numpy.dot(m_lng,m.transpose()) + + mp_lng = mp_lng.reshape(numpy.size(mp_lng),1) + mprime = mprime.reshape(numpy.size(mprime),1) + matrix_mp_lng = numpy.dot(mp_lng,mprime.transpose()) + + f_lng = f_lng.reshape(numpy.size(f_lng),1) + f = f.reshape(numpy.size(f),1) + matrix_f_lng = numpy.dot(f_lng,f.transpose()) + + om_lng = om_lng.reshape(numpy.size(om_lng),1) + omega = omega.reshape(numpy.size(omega),1) + matrix_om_lng = numpy.dot(om_lng,omega.transpose()) + + arg = matrix_d_lng + matrix_m_lng + matrix_mp_lng + matrix_f_lng + matrix_om_lng + + sarg = numpy.sin(arg) + carg = numpy.cos(arg) + + for ii in numpy.arange(n): + nut_long[ii] = 0.0001*numpy.sum((sdelt*t[ii] + sin_lng)*sarg[:,ii]) + nut_obliq[ii] = 0.0001*numpy.sum((cdelt*t[ii] + cos_lng)*carg[:,ii]) + + if numpy.size(jd)==1: + nut_long = nut_long[0] + nut_obliq = nut_obliq[0] + + return nut_long, nut_obliq + + def co_aberration(self,jd,ra,dec): + """ + co_aberration calculates changes to Ra and Dec due to "the effect of aberration". + + Parameters + ---------- + jd = Julian Date (Scalar or vector). + ra = A scalar o vector giving the Right Ascention of interest. + dec = A scalar o vector giving the Declination of interest. + + Return + ------ + d_ra = The correction to right ascension due to aberration (must be added to ra to + get the correct value). + d_dec = The correction to declination due to aberration (must be added to the dec + to get the correct value). + eps = True obliquity of the ecliptic (in radians). + + Examples + -------- + >> Julian = 2462088.7 + >> Ra = 41.547213 + >> Dec = 49.348483 + >> [d_ra,d_dec,eps] = co_aberration(julian,Ra,Dec) + >> print d_ra, d_dec, eps + [ 30.04441796] [ 6.69837858] [ 0.40904059] + + Modification history + -------------------- + Written by Chris O'Dell , Univ. of Wisconsin, June 2002. + Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009. + """ + + # Julian centuries from J2000 of jd. + T = (jd - 2451545.0)/36525.0 + + # Getting obliquity of ecliptic + njd = numpy.size(jd) + jd = numpy.atleast_1d(jd) + ra = numpy.atleast_1d(ra) + dec = numpy.atleast_1d(dec) + + d_psi = numpy.zeros(njd) + d_epsilon = d_psi + for ii in numpy.arange(njd): + [dp,de] = self.nutate(jd[ii]) + d_psi[ii] = dp + d_epsilon[ii] = de + + coeff = 23 + 26/60. + 21.488/3600. + eps0 = coeff*3600. - 46.8150*T - 0.00059*T**2. + 0.001813*T**3. + # True obliquity of the ecliptic in radians + eps = (eps0 + d_epsilon)/3600*Misc_Routines.CoFactors.d2r + + celestialbodies = CelestialBodies() + [sunra,sundec,sunlon,sunobliq] = celestialbodies.sunpos(jd) + + # Earth's orbital eccentricity + e = 0.016708634 - 0.000042037*T - 0.0000001267*T**2. + + # longitude of perihelion, in degrees + pi = 102.93735 + 1.71946*T + 0.00046*T**2 + + # Constant of aberration, in arcseconds + k = 20.49552 + + cd = numpy.cos(dec*Misc_Routines.CoFactors.d2r) ; sd = numpy.sin(dec*Misc_Routines.CoFactors.d2r) + ce = numpy.cos(eps) ; te = numpy.tan(eps) + cp = numpy.cos(pi*Misc_Routines.CoFactors.d2r) ; sp = numpy.sin(pi*Misc_Routines.CoFactors.d2r) + cs = numpy.cos(sunlon*Misc_Routines.CoFactors.d2r) ; ss = numpy.sin(sunlon*Misc_Routines.CoFactors.d2r) + ca = numpy.cos(ra*Misc_Routines.CoFactors.d2r) ; sa = numpy.sin(ra*Misc_Routines.CoFactors.d2r) + + term1 = (ca*cs*ce + sa*ss)/cd + term2 = (ca*cp*ce + sa*sp)/cd + term3 = (cs*ce*(te*cd - sa*sd) + ca*sd*ss) + term4 = (cp*ce*(te*cd - sa*sd) + ca*sd*sp) + + d_ra = -k*term1 + e*k*term2 + d_dec = -k*term3 + e*k*term4 + + return d_ra, d_dec, eps + + def precess(self,ra,dec,equinox1=None,equinox2=None,FK4=0,rad=0): + """ + precess coordinates from EQUINOX1 to EQUINOX2 + + Parameters + ----------- + ra = A scalar o vector giving the Right Ascention of interest. + dec = A scalar o vector giving the Declination of interest. + equinox1 = Original equinox of coordinates, numeric scalar. If omitted, the __Pre- + cess will query for equinox1 and equinox2. + equinox2 = Original equinox of coordinates. + FK4 = If this keyword is set and non-zero, the FK4 (B1950) system will be used + otherwise FK5 (J2000) will be used instead. + rad = If this keyword is set and non-zero, then the input and output RAD and DEC + vectors are in radian rather than degree. + + Return + ------ + ra = Right ascension after precession (scalar or vector) in degrees, unless the rad + keyword is set. + dec = Declination after precession (scalar or vector) in degrees, unless the rad + keyword is set. + + Examples + -------- + >> Ra = 329.88772 + >> Dec = -56.992515 + >> [p_ra,p_dec] = precess(Ra,Dec,1950,1975,FK4=1) + >> print p_ra, p_dec + [ 330.31442971] [-56.87186154] + + Modification history + -------------------- + Written by Wayne Landsman, STI Corporation, August 1986. + Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009. + """ + + npts = numpy.size(ra) + ra = numpy.atleast_1d(ra) + dec = numpy.atleast_1d(dec) + + if rad==0: + ra_rad = ra*Misc_Routines.CoFactors.d2r + dec_rad = dec*Misc_Routines.CoFactors.d2r + else: + ra_rad = ra + dec_rad = dec + + x = numpy.zeros((npts,3)) + x[:,0] = numpy.cos(dec_rad)*numpy.cos(ra_rad) + x[:,1] = numpy.cos(dec_rad)*numpy.sin(ra_rad) + x[:,2] = numpy.sin(dec_rad) + + # Use premat function to get precession matrix from equinox1 to equinox2 + r = self.premat(equinox1,equinox2,FK4) + + x2 = numpy.dot(r,x.transpose()) + + ra_rad = numpy.arctan2(x2[1,:],x2[0,:]) + dec_rad = numpy.arcsin(x2[2,:]) + + if rad==0: + ra = ra_rad/Misc_Routines.CoFactors.d2r + ra = ra + (ra<0)*360. + dec = dec_rad/Misc_Routines.CoFactors.d2r + else: + ra = ra_rad + ra = ra + (ra<0)*numpy.pi*2. + dec = dec_rad + + return ra, dec + + def premat(self,equinox1,equinox2,FK4=0): + """ + premat returns the precession matrix needed to go from EQUINOX1 to EQUINOX2. + + Parameters + ---------- + equinox1 = Original equinox of coordinates, numeric scalar. + equinox2 = Equinox of precessed coordinates. + FK4 = If this keyword is set and non-zero, the FK4 (B1950) system precession angles + are used to compute the precession matrix. The default is to use FK5 (J2000) pre- + cession angles. + + Return + ------ + r = Precession matrix, used to precess equatorial rectangular coordinates. + + Examples + -------- + >> matrix = premat(1950.0,1975.0,FK4=1) + >> print matrix + [[ 9.99981438e-01 -5.58774959e-03 -2.42908517e-03] + [ 5.58774959e-03 9.99984388e-01 -6.78691471e-06] + [ 2.42908517e-03 -6.78633095e-06 9.99997050e-01]] + + Modification history + -------------------- + Written by Wayne Landsman, HSTX Corporation, June 1994. + Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009. + """ + + t = 0.001*(equinox2 - equinox1) + + if FK4==0: + st=0.001*(equinox1 - 2000.) + # Computing 3 rotation angles. + A=Misc_Routines.CoFactors.s2r*t*(23062.181+st*(139.656+0.0139*st)+t*(30.188-0.344*st+17.998*t)) + B=Misc_Routines.CoFactors.s2r*t*t*(79.280+0.410*st+0.205*t)+A + C=Misc_Routines.CoFactors.s2r*t*(20043.109-st*(85.33+0.217*st)+ t*(-42.665-0.217*st-41.833*t)) + else: + st=0.001*(equinox1 - 1900) + # Computing 3 rotation angles + A=Misc_Routines.CoFactors.s2r*t*(23042.53+st*(139.75+0.06*st)+t*(30.23-0.27*st+18.0*t)) + B=Misc_Routines.CoFactors.s2r*t*t*(79.27+0.66*st+0.32*t)+A + C=Misc_Routines.CoFactors.s2r*t*(20046.85-st*(85.33+0.37*st)+t*(-42.67-0.37*st-41.8*t)) + + sina = numpy.sin(A); sinb = numpy.sin(B); sinc = numpy.sin(C) + cosa = numpy.cos(A); cosb = numpy.cos(B); cosc = numpy.cos(C) + + r = numpy.zeros((3,3)) + r[:,0] = numpy.array([cosa*cosb*cosc-sina*sinb,sina*cosb+cosa*sinb*cosc,cosa*sinc]) + r[:,1] = numpy.array([-cosa*sinb-sina*cosb*cosc,cosa*cosb-sina*sinb*cosc,-sina*sinc]) + r[:,2] = numpy.array([-cosb*sinc,-sinb*sinc,cosc]) + + return r + + def cirrange(self,angle,rad=0): + """ + cirrange forces an angle into the range 0<= angle < 360. + + Parameters + ---------- + angle = The angle to modify, in degrees. Can be scalar or vector. + rad = Set to 1 if the angle is specified in radians rather than degrees. It is for- + ced into the range 0 <= angle < 2 PI + + Return + ------ + angle = The angle after the modification. + + Example + ------- + >> angle = cirrange(numpy.array([420,400,361])) + >> print angle + >> [60, 40, 1] + + Modification History + -------------------- + Written by Michael R. Greason, Hughes STX, 10 February 1994. + Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009. + """ + + angle = numpy.atleast_1d(angle) + + if rad==1: + cnst = numpy.pi*2. + elif rad==0: + cnst = 360. + + # Deal with the lower limit. + angle = angle % cnst + + # Deal with negative values, if way + neg = numpy.where(angle<0.0) + if neg[0].size>0: angle[neg] = angle[neg] + cnst + + return angle + + +class CelestialBodies(EquatorialCorrections): + def __init__(self): + """ + CelestialBodies class creates a object to call methods of celestial bodies location. + + Modification History + -------------------- + Converted to Object-oriented Programming by Freddy Galindo, ROJ, 27 September 2009. + """ + + EquatorialCorrections.__init__(self) + + def sunpos(self,jd,rad=0): + """ + sunpos method computes the RA and Dec of the Sun at a given date. + + Parameters + ---------- + jd = The julian date of the day (and time), scalar or vector. + rad = If this keyword is set and non-zero, then the input and output RAD and DEC + vectors are in radian rather than degree. + + Return + ------ + ra = The right ascension of the sun at that date in degrees. + dec = The declination of the sun at that date in degrees. + elong = Ecliptic longitude of the sun at that date in degrees. + obliquity = The declination of the sun at that date in degrees. + + Examples + -------- + >> jd = 2466880 + >> [ra,dec,elong,obliquity] = sunpos(jd) + >> print ra, dec, elong, obliquity + [ 275.53499556] [-23.33840558] [ 275.08917968] [ 23.43596165] + + >> [ra,dec,elong,obliquity] = sunpos(jd,rad=1) + >> print ra, dec, elong, obliquity + [ 4.80899288] [-0.40733202] [ 4.80121192] [ 0.40903469] + + >> jd = 2450449.5 + numpy.arange(365) + >> [ra,dec,elong,obliquity] = sunpos(jd) + + Modification history + -------------------- + Written by Micheal R. Greason, STX Corporation, 28 October 1988. + Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009. + """ + + jd = numpy.atleast_1d(jd) + + # Form time in Julian centuries from 1900. + t = (jd -2415020.0)/36525.0 + + # Form sun's mean longitude + l = (279.696678+((36000.768925*t) % 360.0))*3600.0 + + # Allow for ellipticity of the orbit (equation of centre) using the Earth's mean + # anomoly ME + me = 358.475844 + ((35999.049750*t) % 360.0) + ellcor = (6910.1 - 17.2*t)*numpy.sin(me*Misc_Routines.CoFactors.d2r) + 72.3*numpy.sin(2.0*me*Misc_Routines.CoFactors.d2r) + l = l + ellcor + + # Allow for the Venus perturbations using the mean anomaly of Venus MV + mv = 212.603219 + ((58517.803875*t) % 360.0) + vencorr = 4.8*numpy.cos((299.1017 + mv - me)*Misc_Routines.CoFactors.d2r) + \ + 5.5*numpy.cos((148.3133 + 2.0*mv - 2.0*me )*Misc_Routines.CoFactors.d2r) + \ + 2.5*numpy.cos((315.9433 + 2.0*mv - 3.0*me )*Misc_Routines.CoFactors.d2r) + \ + 1.6*numpy.cos((345.2533 + 3.0*mv - 4.0*me )*Misc_Routines.CoFactors.d2r) + \ + 1.0*numpy.cos((318.15 + 3.0*mv - 5.0*me )*Misc_Routines.CoFactors.d2r) + l = l + vencorr + + # Allow for the Mars perturbations using the mean anomaly of Mars MM + mm = 319.529425 + ((19139.858500*t) % 360.0) + marscorr = 2.0*numpy.cos((343.8883 - 2.0*mm + 2.0*me)*Misc_Routines.CoFactors.d2r ) + \ + 1.8*numpy.cos((200.4017 - 2.0*mm + me)*Misc_Routines.CoFactors.d2r) + l = l + marscorr + + # Allow for the Jupiter perturbations using the mean anomaly of Jupiter MJ + mj = 225.328328 + ((3034.6920239*t) % 360.0) + jupcorr = 7.2*numpy.cos((179.5317 - mj + me )*Misc_Routines.CoFactors.d2r) + \ + 2.6*numpy.cos((263.2167 - mj)*Misc_Routines.CoFactors.d2r) + \ + 2.7*numpy.cos((87.1450 - 2.0*mj + 2.0*me)*Misc_Routines.CoFactors.d2r) + \ + 1.6*numpy.cos((109.4933 - 2.0*mj + me)*Misc_Routines.CoFactors.d2r) + l = l + jupcorr + + # Allow for Moons perturbations using mean elongation of the Moon from the Sun D + d = 350.7376814 + ((445267.11422*t) % 360.0) + mooncorr = 6.5*numpy.sin(d*Misc_Routines.CoFactors.d2r) + l = l + mooncorr + + # Allow for long period terms + longterm = + 6.4*numpy.sin((231.19 + 20.20*t)*Misc_Routines.CoFactors.d2r) + l = l + longterm + l = (l + 2592000.0) % 1296000.0 + longmed = l/3600.0 + + # Allow for Aberration + l = l - 20.5 + + # Allow for Nutation using the longitude of the Moons mean node OMEGA + omega = 259.183275 - ((1934.142008*t) % 360.0) + l = l - 17.2*numpy.sin(omega*Misc_Routines.CoFactors.d2r) + + # Form the True Obliquity + oblt = 23.452294 - 0.0130125*t + (9.2*numpy.cos(omega*Misc_Routines.CoFactors.d2r))/3600.0 + + # Form Right Ascension and Declination + l = l/3600.0 + ra = numpy.arctan2((numpy.sin(l*Misc_Routines.CoFactors.d2r)*numpy.cos(oblt*Misc_Routines.CoFactors.d2r)),numpy.cos(l*Misc_Routines.CoFactors.d2r)) + + neg = numpy.where(ra < 0.0) + if neg[0].size > 0: ra[neg] = ra[neg] + 2.0*numpy.pi + + dec = numpy.arcsin(numpy.sin(l*Misc_Routines.CoFactors.d2r)*numpy.sin(oblt*Misc_Routines.CoFactors.d2r)) + + if rad==1: + oblt = oblt*Misc_Routines.CoFactors.d2r + longmed = longmed*Misc_Routines.CoFactors.d2r + else: + ra = ra/Misc_Routines.CoFactors.d2r + dec = dec/Misc_Routines.CoFactors.d2r + + return ra, dec, longmed, oblt + + def moonpos(self,jd,rad=0): + """ + moonpos method computes the RA and Dec of the Moon at specified Julian date(s). + + Parameters + ---------- + jd = The julian date of the day (and time), scalar or vector. + rad = If this keyword is set and non-zero, then the input and output RAD and DEC + vectors are in radian rather than degree. + + Return + ------ + ra = The right ascension of the sun at that date in degrees. + dec = The declination of the sun at that date in degrees. + dist = The Earth-moon distance in kilometers (between the center of the Earth and + the center of the moon). + geolon = Apparent longitude of the moon in degrees, referred to the ecliptic of the + specified date(s). + geolat = Apparent latitude the moon in degrees, referred to the ecliptic of the + specified date(s). + + Examples + -------- + >> jd = 2448724.5 + >> [ra,dec,dist,geolon,geolat] = sunpos(jd) + >> print ra, dec, dist, geolon, geolat + [ 134.68846855] [ 13.76836663] [ 368409.68481613] [ 133.16726428] [-3.22912642] + + >> [ra,dec,dist,geolon, geolat] = sunpos(jd,rad=1) + >> print ra, dec, dist, geolon, geolat + [ 2.35075724] [ 0.24030333] [ 368409.68481613] [ 2.32420722] [-0.05635889] + + >> jd = 2450449.5 + numpy.arange(365) + >> [ra,dec,dist,geolon, geolat] = sunpos(jd) + + Modification history + -------------------- + Written by Micheal R. Greason, STX Corporation, 31 October 1988. + Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. + """ + + jd = numpy.atleast_1d(jd) + + # Form time in Julian centuries from 1900. + t = (jd - 2451545.0)/36525.0 + + d_lng = numpy.array([0,2,2,0,0,0,2,2,2,2,0,1,0,2,0,0,4,0,4,2,2,1,1,2,2,4,2,0,2,2,1,2,\ + 0,0,2,2,2,4,0,3,2,4,0,2,2,2,4,0,4,1,2,0,1,3,4,2,0,1,2,2]) + + m_lng = numpy.array([0,0,0,0,1,0,0,-1,0,-1,1,0,1,0,0,0,0,0,0,1,1,0,1,-1,0,0,0,1,0,-1,\ + 0,-2,1,2,-2,0,0,-1,0,0,1,-1,2,2,1,-1,0,0,-1,0,1,0,1,0,0,-1,2,1,0,0]) + + mp_lng = numpy.array([1,-1,0,2,0,0,-2,-1,1,0,-1,0,1,0,1,1,-1,3,-2,-1,0,-1,0,1,2,0,-3,\ + -2,-1,-2,1,0,2,0,-1,1,0,-1,2,-1,1,-2,-1,-1,-2,0,1,4,0,-2,0,2,1,-2,-3,2,1,-1,3,-1]) + + f_lng = numpy.array([0,0,0,0,0,2,0,0,0,0,0,0,0,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,\ + 0,0,0,0,-2,2,0,2,0,0,0,0,0,0,-2,0,0,0,0,-2,-2,0,0,0,0,0,0,0,-2]) + + sin_lng = numpy.array([6288774,1274027,658314,213618,-185116,-114332,58793,57066,\ + 53322,45758,-40923,-34720,-30383,15327,-12528,10980,10675,10034,8548,-7888,\ + -6766,-5163,4987,4036,3994,3861,3665,-2689,-2602,2390,-2348,2236,-2120,-2069,\ + 2048,-1773,-1595,1215,-1110,-892,-810,759,-713,-700,691,596,549,537,520,-487,\ + -399,-381,351,-340,330,327,-323,299,294,0.0]) + + cos_lng = numpy.array([-20905355,-3699111,-2955968,-569925,48888,-3149,246158,-152138,\ + -170733,-204586,-129620,108743,104755,10321,0,79661,-34782,-23210,-21636,24208,\ + 30824,-8379,-16675,-12831,-10445,-11650,14403,-7003,0,10056,6322, -9884,5751,0,\ + -4950,4130,0,-3958,0,3258,2616,-1897,-2117,2354,0,0,-1423,-1117,-1571,-1739,0, \ + -4421,0,0,0,0,1165,0,0,8752.0]) + + d_lat = numpy.array([0,0,0,2,2,2,2,0,2,0,2,2,2,2,2,2,2,0,4,0,0,0,1,0,0,0,1,0,4,4,0,4,\ + 2,2,2,2,0,2,2,2,2,4,2,2,0,2,1,1,0,2,1,2,0,4,4,1,4,1,4,2]) + + m_lat = numpy.array([0,0,0,0,0,0,0,0,0,0,-1,0,0,1,-1,-1,-1,1,0,1,0,1,0,1,1,1,0,0,0,0,\ + 0,0,0,0,-1,0,0,0,0,1,1,0,-1,-2,0,1,1,1,1,1,0,-1,1,0,-1,0,0,0,-1,-2]) + + mp_lat = numpy.array([0,1,1,0,-1,-1,0,2,1,2,0,-2,1,0,-1,0,-1,-1,-1,0,0,-1,0,1,1,0,0,\ + 3,0,-1,1,-2,0,2,1,-2,3,2,-3,-1,0,0,1,0,1,1,0,0,-2,-1,1,-2,2,-2,-1,1,1,-1,0,0]) + + f_lat = numpy.array([1,1,-1,-1,1,-1,1,1,-1,-1,-1,-1,1,-1,1,1,-1,-1,-1,1,3,1,1,1,-1,\ + -1,-1,1,-1,1,-3,1,-3,-1,-1,1,-1,1,-1,1,1,1,1,-1,3,-1,-1,1,-1,-1,1,-1,1,-1,-1, \ + -1,-1,-1,-1,1]) + + sin_lat = numpy.array([5128122,280602,277693,173237,55413,46271, 32573, 17198, 9266, \ + 8822,8216,4324,4200,-3359,2463,2211,2065,-1870,1828,-1794, -1749, -1565, -1491, \ + -1475,-1410,-1344,-1335,1107,1021,833,777,671,607,596,491,-451,439,422,421,-366,\ + -351,331,315,302,-283,-229,223,223,-220,-220,-185,181,-177,176, 166, -164, 132, \ + -119,115,107.0]) + + # Mean longitude of the moon refered to mean equinox of the date. + coeff0 = numpy.array([-1./6.5194e7,1./538841.,-0.0015786,481267.88123421,218.3164477]) + lprimed = numpy.poly1d(coeff0) + lprimed = lprimed(t) + lprimed = self.cirrange(lprimed,rad=0) + lprime = lprimed*Misc_Routines.CoFactors.d2r + + # Mean elongation of the moon + coeff1 = numpy.array([-1./1.13065e8,1./545868.,-0.0018819,445267.1114034,297.8501921]) + d = numpy.poly1d(coeff1) + d = d(t)*Misc_Routines.CoFactors.d2r + d = self.cirrange(d,rad=1) + + # Sun's mean anomaly + coeff2 = numpy.array([1.0/2.449e7,-0.0001536,35999.0502909,357.5291092]) + M = numpy.poly1d(coeff2) + M = M(t)*Misc_Routines.CoFactors.d2r + M = self.cirrange(M,rad=1) + + # Moon's mean anomaly + coeff3 = numpy.array([-1.0/1.4712e7,1.0/6.9699e4,0.0087414,477198.8675055,134.9633964]) + Mprime = numpy.poly1d(coeff3) + Mprime = Mprime(t)*Misc_Routines.CoFactors.d2r + Mprime = self.cirrange(Mprime,rad=1) + + # Moon's argument of latitude + coeff4 = numpy.array([1.0/8.6331e8,-1.0/3.526e7,-0.0036539,483202.0175233,93.2720950]) + F = numpy.poly1d(coeff4) + F = F(t)*Misc_Routines.CoFactors.d2r + F = self.cirrange(F,rad=1) + + # Eccentricity of Earth's orbit around the sun + e = 1 - 0.002516*t - 7.4e-6*(t**2.) + e2 = e**2. + + ecorr1 = numpy.where((numpy.abs(m_lng))==1) + ecorr2 = numpy.where((numpy.abs(m_lat))==1) + ecorr3 = numpy.where((numpy.abs(m_lng))==2) + ecorr4 = numpy.where((numpy.abs(m_lat))==2) + + # Additional arguments. + A1 = (119.75 + 131.849*t)*Misc_Routines.CoFactors.d2r + A2 = (53.09 + 479264.290*t)*Misc_Routines.CoFactors.d2r + A3 = (313.45 + 481266.484*t)*Misc_Routines.CoFactors.d2r + suml_add = 3958.*numpy.sin(A1) + 1962.*numpy.sin(lprime - F) + 318*numpy.sin(A2) + sumb_add = -2235.*numpy.sin(lprime) + 382.*numpy.sin(A3) + 175.*numpy.sin(A1-F) + \ + 175.*numpy.sin(A1 + F) + 127.*numpy.sin(lprime - Mprime) - 115.*numpy.sin(lprime + Mprime) + + # Sum the periodic terms + geolon = numpy.zeros(jd.size) + geolat = numpy.zeros(jd.size) + dist = numpy.zeros(jd.size) + + for i in numpy.arange(jd.size): + sinlng = sin_lng + coslng = cos_lng + sinlat = sin_lat + + sinlng[ecorr1] = e[i]*sinlng[ecorr1] + coslng[ecorr1] = e[i]*coslng[ecorr1] + sinlat[ecorr2] = e[i]*sinlat[ecorr2] + sinlng[ecorr3] = e2[i]*sinlng[ecorr3] + coslng[ecorr3] = e2[i]*coslng[ecorr3] + sinlat[ecorr4] = e2[i]*sinlat[ecorr4] + + arg = d_lng*d[i] + m_lng*M[i] + mp_lng*Mprime[i] + f_lng*F[i] + geolon[i] = lprimed[i] + (numpy.sum(sinlng*numpy.sin(arg)) + suml_add[i] )/1.e6 + dist[i] = 385000.56 + numpy.sum(coslng*numpy.cos(arg))/1.e3 + arg = d_lat*d[i] + m_lat*M[i] + mp_lat*Mprime[i] + f_lat*F[i] + geolat[i] = (numpy.sum(sinlat*numpy.sin(arg)) + sumb_add[i])/1.e6 + + [nlon, elon] = self.nutate(jd) + geolon = geolon + nlon/3.6e3 + geolon = self.cirrange(geolon,rad=0) + lamb = geolon*Misc_Routines.CoFactors.d2r + beta = geolat*Misc_Routines.CoFactors.d2r + + # Find mean obliquity and convert lamb, beta to RA, Dec + c = numpy.array([2.45,5.79,27.87,7.12,-39.05,-249.67,-51.38,1999.25,-1.55,-4680.93, \ + 21.448]) + junk = numpy.poly1d(c); + epsilon = 23. + (26./60.) + (junk(t/1.e2)/3600.) + # True obliquity in radians + eps = (epsilon + elon/3600. )*Misc_Routines.CoFactors.d2r + + ra = numpy.arctan2(numpy.sin(lamb)*numpy.cos(eps)-numpy.tan(beta)*numpy.sin(eps),numpy.cos(lamb)) + ra = self.cirrange(ra,rad=1) + + dec = numpy.arcsin(numpy.sin(beta)*numpy.cos(eps) + numpy.cos(beta)*numpy.sin(eps)*numpy.sin(lamb)) + + if rad==1: + geolon = lamb + geolat = beta + else: + ra = ra/Misc_Routines.CoFactors.d2r + dec = dec/Misc_Routines.CoFactors.d2r + + return ra, dec, dist, geolon, geolat + + def hydrapos(self): + """ + hydrapos method returns RA and Dec provided by Bill Coles (Oct 2003). + + Parameters + ---------- + None + + Return + ------ + ra = The right ascension of the sun at that date in degrees. + dec = The declination of the sun at that date in degrees. + Examples + -------- + >> [ra,dec] = hydrapos() + >> print ra, dec + 139.45 -12.0833333333 + + Modification history + -------------------- + Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. + """ + + ra = (9. + 17.8/60.)*15. + dec = -(12. + 5./60.) + + return ra, dec + + + def skynoise_jro(self,dec_cut=-11.95,filename='skynoise_jro.dat',filepath=None): + """ + hydrapos returns RA and Dec provided by Bill Coles (Oct 2003). + + Parameters + ---------- + dec_cut = A scalar giving the declination to get a cut of the skynoise over Jica- + marca. The default value is -11.95. + filename = A string to specify name the skynoise file. The default value is skynoi- + se_jro.dat + + Return + ------ + maxra = The maximum right ascension to the declination used to get a cut. + ra = The right ascension. + Examples + -------- + >> [maxra,ra] = skynoise_jro() + >> print maxra, ra + 139.45 -12.0833333333 + + Modification history + -------------------- + Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. + """ + + if filepath==None: + filepath = '/app/utils/' + + f = open(os.path.join(filepath,filename),'rb') + + # Reading SkyNoise Power (lineal scale) + ha_sky = numpy.fromfile(f,numpy.dtype([('var','> [maxra,ra] = skynoise_jro() + >> print maxra, ra + 139.45 -12.0833333333 + + Modification history + -------------------- + Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. + """ + + # Defining date to compute SkyNoise. + [year, month, dom, hour, mis, secs] = TimeTools.Julian(jd).change2time() + is_dom = (month==9) & (dom==21) + if is_dom: + tmp = jd + jd = TimeTools.Time(year,9,22).change2julian() + dom = 22 + + # Reading SkyNoise + if filepath==None:filepath='./resource' + f = open(os.path.join(filepath,filename)) + + lines = f.read() + f.close() + + nlines = 99 + lines = lines.split('\n') + data = numpy.zeros((2,nlines))*numpy.float32(0.) + for ii in numpy.arange(nlines): + line = numpy.array([lines[ii][0:6],lines[ii][6:]]) + data[:,ii] = numpy.float32(line) + + # Getting SkyNoise to the date desired. + otime = data[0,:]*60.0 + opowr = data[1,:] + + hour = numpy.array([0,23]); + mins = numpy.array([0,59]); + secs = numpy.array([0,59]); + LTrange = TimeTools.Time(year,month,dom,hour,mins,secs).change2julday() + LTtime = LTrange[0] + numpy.arange(1440)*((LTrange[1] - LTrange[0])/(1440.-1)) + lst = TimeTools.Julian(LTtime + (-3600.*ut/86400.)).change2lst() + + ipowr = lst*0.0 + # Interpolating using scipy (inside max and min "x") + otime = otime/3600. + val = numpy.where((lst>numpy.min(otime)) & (lstnumpy.max(otime)) + if uval[0].size>0: + ii = numpy.min(uval[0]) + m = (ipowr[ii-1] - ipowr[ii-2])/(lst[ii-1] - lst[ii-2]) + b = ipowr[ii-1] - m*lst[ii-1] + ipowr[uval] = m*lst[uval] + b + + if is_dom: + lst = numpy.roll(lst,4) + ipowr = numpy.roll(ipowr,4) + + new_lst = numpy.int32(lst*3600.) + new_pow = ipowr + + return ipowr, LTtime, lst + + +class AltAz(EquatorialCorrections): + def __init__(self,alt,az,jd,lat=-11.95,lon=-76.8667,WS=0,altitude=500,nutate_=0,precess_=0,\ + aberration_=0,B1950=0): + """ + The AltAz class creates an object which represents the target position in horizontal + coordinates (alt-az) and allows to convert (using the methods) from this coordinate + system to others (e.g. Equatorial). + + Parameters + ---------- + alt = Altitude in degrees. Scalar or vector. + az = Azimuth angle in degrees (measured EAST from NORTH, but see keyword WS). Sca- + lar or vector. + jd = Julian date. Scalar or vector. + lat = North geodetic latitude of location in degrees. The default value is -11.95. + lon = East longitude of location in degrees. The default value is -76.8667. + WS = Set this to 1 to get the azimuth measured westward from south. + altitude = The altitude of the observing location, in meters. The default 500. + nutate_ = Set this to 1 to force nutation, 0 for no nutation. + precess_ = Set this to 1 to force precession, 0 for no precession. + aberration_ = Set this to 1 to force aberration correction, 0 for no correction. + B1950 = Set this if your RA and DEC are specified in B1950, FK4 coordinates (ins- + tead of J2000, FK5) + + Modification History + -------------------- + Converted to Object-oriented Programming by Freddy Galindo, ROJ, 26 September 2009. + """ + + EquatorialCorrections.__init__(self) + + self.alt = numpy.atleast_1d(alt) + self.az = numpy.atleast_1d(az) + self.jd = numpy.atleast_1d(jd) + self.lat = lat + self.lon = lon + self.WS = WS + self.altitude = altitude + + self.nutate_ = nutate_ + self.aberration_ = aberration_ + self.precess_ = precess_ + self.B1950 = B1950 + + def change2equatorial(self): + """ + change2equatorial method converts horizon (Alt-Az) coordinates to equatorial coordi- + nates (ra-dec). + + Return + ------ + ra = Right ascension of object (J2000) in degrees (FK5). Scalar or vector. + dec = Declination of object (J2000), in degrees (FK5). Scalar or vector. + ha = Hour angle in degrees. + + Example + ------- + >> alt = 88.5401 + >> az = -128.990 + >> jd = 2452640.5 + >> ObjAltAz = AltAz(alt,az,jd) + >> [ra, dec, ha] = ObjAltAz.change2equatorial() + >> print ra, dec, ha + [ 22.20280632] [-12.86610025] [ 1.1638927] + + Modification History + -------------------- + Written Chris O'Dell Univ. of Wisconsin-Madison, May 2002. + Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009. + """ + + az = self.az + alt = self.alt + if self.WS>0:az = az -180. + ra_tmp = numpy.zeros(numpy.size(self.jd)) + 45. + dec_tmp = numpy.zeros(numpy.size(self.jd)) + 45. + [dra1,ddec1,eps,d_psi,d_eps] = self.co_nutate(self.jd,ra_tmp, dec_tmp) + + # Getting local mean sidereal time (lmst) + lmst = TimeTools.Julian(self.jd[0]).change2lst() + lmst = lmst*Misc_Routines.CoFactors.h2d + # Getting local apparent sidereal time (last) + last = lmst + d_psi*numpy.cos(eps)/3600. + + # Now do the spherical trig to get APPARENT hour angle and declination (Degrees). + [ha, dec] = self.change2HaDec() + + # Finding Right Ascension (in degrees, from 0 to 360.) + ra = (last - ha + 360.) % 360. + + # Calculate NUTATION and ABERRATION Correction to Ra-Dec + [dra1, ddec1,eps,d_psi,d_eps] = self.co_nutate(self.jd,ra,dec) + [dra2,ddec2,eps] = self.co_aberration(self.jd,ra,dec) + + # Make Nutation and Aberration correction (if wanted) + ra = ra - (dra1*self.nutate_ + dra2*self.aberration_)/3600. + dec = dec - (ddec1*self.nutate_ + ddec2*self.aberration_)/3600. + + # Computing current equinox + j_now = (self.jd - 2451545.)/365.25 + 2000 + + # Precess coordinates to current date + if self.precess_==1: + njd = numpy.size(self.jd) + for ii in numpy.arange(njd): + ra_i = ra[ii] + dec_i = dec[ii] + now = j_now[ii] + + if self.B1950==1: + [ra_i,dec_i] = self.precess(ra_i,dec_i,now,1950.,FK4=1) + elif self.B1950==0: + [ra_i,dec_i] = self.precess(ra_i,dec_i,now,2000.,FK4=0) + + ra[ii] = ra_i + dec[ii] = dec_i + + return ra, dec, ha + + def change2HaDec(self): + """ + change2HaDec method converts from horizon (Alt-Az) coordinates to hour angle and de- + clination. + + Return + ------ + ha = The local apparent hour angle, in degrees. The hour angle is the time that ri- + ght ascension of 0 hours crosses the local meridian. It is unambiguisoly defined. + dec = The local apparent declination, in degrees. + + Example + ------- + >> alt = 88.5401 + >> az = -128.990 + >> jd = 2452640.5 + >> ObjAltAz = AltAz(alt,az,jd) + >> [ha, dec] = ObjAltAz.change2HaDec() + >> print ha, dec + [ 1.1638927] [-12.86610025] + + Modification History + -------------------- + Written Chris O'Dell Univ. of Wisconsin-Madison, May 2002. + Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009. + """ + + alt_r = numpy.atleast_1d(self.alt*Misc_Routines.CoFactors.d2r) + az_r = numpy.atleast_1d(self.az*Misc_Routines.CoFactors.d2r) + lat_r = numpy.atleast_1d(self.lat*Misc_Routines.CoFactors.d2r) + + # Find local hour angle (in degrees, from 0 to 360.) + y_ha = -1*numpy.sin(az_r)*numpy.cos(alt_r) + x_ha = -1*numpy.cos(az_r)*numpy.sin(lat_r)*numpy.cos(alt_r) + numpy.sin(alt_r)*numpy.cos(lat_r) + + ha = numpy.arctan2(y_ha,x_ha) + ha = ha/Misc_Routines.CoFactors.d2r + + w = numpy.where(ha<0.) + if w[0].size>0:ha[w] = ha[w] + 360. + ha = ha % 360. + + # Find declination (positive if north of celestial equatorial, negative if south) + sindec = numpy.sin(lat_r)*numpy.sin(alt_r) + numpy.cos(lat_r)*numpy.cos(alt_r)*numpy.cos(az_r) + dec = numpy.arcsin(sindec)/Misc_Routines.CoFactors.d2r + + return ha, dec + + +class Equatorial(EquatorialCorrections): + def __init__(self,ra,dec,jd,lat=-11.95,lon=-76.8667,WS=0,altitude=500,nutate_=0,precess_=0,\ + aberration_=0,B1950=0): + """ + The Equatorial class creates an object which represents the target position in equa- + torial coordinates (ha-dec) and allows to convert (using the class methods) from + this coordinate system to others (e.g. AltAz). + + Parameters + ---------- + ra = Right ascension of object (J2000) in degrees (FK5). Scalar or vector. + dec = Declination of object (J2000), in degrees (FK5). Scalar or vector. + jd = Julian date. Scalar or vector. + lat = North geodetic latitude of location in degrees. The default value is -11.95. + lon = East longitude of location in degrees. The default value is -76.8667. + WS = Set this to 1 to get the azimuth measured westward from south. + altitude = The altitude of the observing location, in meters. The default 500. + nutate = Set this to 1 to force nutation, 0 for no nutation. + precess = Set this to 1 to force precession, 0 for no precession. + aberration = Set this to 1 to force aberration correction, 0 for no correction. + B1950 = Set this if your RA and DEC are specified in B1950, FK4 coordinates (ins- + tead of J2000, FK5) + + Modification History + -------------------- + Converted to Object-oriented Programming by Freddy Galindo, ROJ, 29 September 2009. + """ + + EquatorialCorrections.__init__(self) + + self.ra = numpy.atleast_1d(ra) + self.dec = numpy.atleast_1d(dec) + self.jd = numpy.atleast_1d(jd) + self.lat = lat + self.lon = lon + self.WS = WS + self.altitude = altitude + + self.nutate_ = nutate_ + self.aberration_ = aberration_ + self.precess_ = precess_ + self.B1950 = B1950 + + def change2AltAz(self): + """ + change2AltAz method converts from equatorial coordinates (ha-dec) to horizon coordi- + nates (alt-az). + + Return + ------ + alt = Altitude in degrees. Scalar or vector. + az = Azimuth angle in degrees (measured EAST from NORTH, but see keyword WS). Sca- + lar or vector. + ha = Hour angle in degrees. + + Example + ------- + >> ra = 43.370609 + >> dec = -28.0000 + >> jd = 2452640.5 + >> ObjEq = Equatorial(ra,dec,jd) + >> [alt, az, ha] = ObjEq.change2AltAz() + >> print alt, az, ha + [ 65.3546497] [ 133.58753124] [ 339.99609002] + + Modification History + -------------------- + Written Chris O'Dell Univ. of Wisconsin-Madison. May 2002 + Converted to Python by Freddy R. Galindo, ROJ, 29 September 2009. + """ + + ra = self.ra + dec = self.dec + + # Computing current equinox + j_now = (self.jd - 2451545.)/365.25 + 2000 + + # Precess coordinates to current date + if self.precess_==1: + njd = numpy.size(self.jd) + for ii in numpy.arange(njd): + ra_i = ra[ii] + dec_i = dec[ii] + now = j_now[ii] + + if self.B1950==1: + [ra_i,dec_i] = self.precess(ra_i,dec_i,now,1950.,FK4=1) + elif self.B1950==0: + [ra_i,dec_i] = self.precess(ra_i,dec_i,now,2000.,FK4=0) + + ra[ii] = ra_i + dec[ii] = dec_i + + # Calculate NUTATION and ABERRATION Correction to Ra-Dec + [dra1, ddec1,eps,d_psi,d_eps] = self.co_nutate(self.jd,ra,dec) + [dra2,ddec2,eps] = self.co_aberration(self.jd,ra,dec) + + # Make Nutation and Aberration correction (if wanted) + ra = ra + (dra1*self.nutate_ + dra2*self.aberration_)/3600. + dec = dec + (ddec1*self.nutate_ + ddec2*self.aberration_)/3600. + + # Getting local mean sidereal time (lmst) + lmst = TimeTools.Julian(self.jd).change2lst() + + lmst = lmst*Misc_Routines.CoFactors.h2d + # Getting local apparent sidereal time (last) + last = lmst + d_psi*numpy.cos(eps)/3600. + + # Finding Hour Angle (in degrees, from 0 to 360.) + ha = last - ra + w = numpy.where(ha<0.) + if w[0].size>0:ha[w] = ha[w] + 360. + ha = ha % 360. + + # Now do the spherical trig to get APPARENT hour angle and declination (Degrees). + [alt, az] = self.HaDec2AltAz(ha,dec) + + return alt, az, ha + + def HaDec2AltAz(self,ha,dec): + """ + HaDec2AltAz convert hour angle and declination (ha-dec) to horizon coords (alt-az). + + Parameters + ---------- + ha = The local apparent hour angle, in DEGREES, scalar or vector. + dec = The local apparent declination, in DEGREES, scalar or vector. + + Return + ------ + alt = Altitude in degrees. Scalar or vector. + az = Azimuth angle in degrees (measured EAST from NORTH, but see keyword WS). Sca- + lar or vector. + + Modification History + -------------------- + Written Chris O'Dell Univ. of Wisconsin-Madison, May 2002. + Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009. + """ + + sh = numpy.sin(ha*Misc_Routines.CoFactors.d2r) ; ch = numpy.cos(ha*Misc_Routines.CoFactors.d2r) + sd = numpy.sin(dec*Misc_Routines.CoFactors.d2r) ; cd = numpy.cos(dec*Misc_Routines.CoFactors.d2r) + sl = numpy.sin(self.lat*Misc_Routines.CoFactors.d2r) ; cl = numpy.cos(self.lat*Misc_Routines.CoFactors.d2r) + + x = -1*ch*cd*sl + sd*cl + y = -1*sh*cd + z = ch*cd*cl + sd*sl + r = numpy.sqrt(x**2. + y**2.) + + az = numpy.arctan2(y,x)/Misc_Routines.CoFactors.d2r + alt = numpy.arctan2(z,r)/Misc_Routines.CoFactors.d2r + + # correct for negative az. + w = numpy.where(az<0.) + if w[0].size>0:az[w] = az[w] + 360. + + # Convert az to West from South, if desired + if self.WS==1: az = (az + 180.) % 360. + + return alt, az + + +class Geodetic(): + def __init__(self,lat=-11.95,alt=0): + """ + The Geodetic class creates an object which represents the real position on earth of + a target (Geodetic Coordinates: lat-alt) and allows to convert (using the class me- + thods) from this coordinate system to others (e.g. geocentric). + + Parameters + ---------- + lat = Geodetic latitude of location in degrees. The default value is -11.95. + + alt = Geodetic altitude (km). The default value is 0. + + Modification History + -------------------- + Converted to Object-oriented Programming by Freddy R. Galindo, ROJ, 02 October 2009. + """ + + self.lat = numpy.atleast_1d(lat) + self.alt = numpy.atleast_1d(alt) + + self.a = 6378.16 + self.ab2 = 1.0067397 + self.ep2 = 0.0067397 + + def change2geocentric(self): + """ + change2geocentric method converts from Geodetic to Geocentric coordinates. The re- + ference geoid is that adopted by the IAU in 1964. + + Return + ------ + gclat = Geocentric latitude (in degrees), scalar or vector. + gcalt = Geocentric radial distance (km), scalar or vector. + + Example + ------- + >> ObjGeoid = Geodetic(lat=-11.95,alt=0) + >> [gclat, gcalt] = ObjGeoid.change2geocentric() + >> print gclat, gcalt + [-11.87227742] [ 6377.25048195] + + Modification History + -------------------- + Converted to Python by Freddy R. Galindo, ROJ, 02 October 2009. + """ + + gdl = self.lat*Misc_Routines.CoFactors.d2r + slat = numpy.sin(gdl) + clat = numpy.cos(gdl) + slat2 = slat**2. + clat2 = (self.ab2*clat)**2. + + sbet = slat/numpy.sqrt(slat2 + clat2) + sbet2 = (sbet**2.) # < 1 + noval = numpy.where(sbet2>1) + if noval[0].size>0:sbet2[noval] = 1 + cbet = numpy.sqrt(1. - sbet2) + + rgeoid = self.a/numpy.sqrt(1. + self.ep2*sbet2) + + x = rgeoid*cbet + self.alt*clat + y = rgeoid*sbet + self.alt*slat + + gcalt = numpy.sqrt(x**2. + y**2.) + gclat = numpy.arctan2(y,x)/Misc_Routines.CoFactors.d2r + + return gclat, gcalt diff --git a/schainpy/model/utils/BField.py b/schainpy/model/utils/BField.py new file mode 100644 index 0000000..7608c9c --- /dev/null +++ b/schainpy/model/utils/BField.py @@ -0,0 +1,507 @@ +""" + +Converted to Object-oriented Programming by Freddy Galindo, ROJ, 07 October 2009. +Added to signal Chain by Joab Apaza, ROJ, Jun 2023. +""" + +import numpy +import time +import os +from scipy.special import lpmn +from schainpy.model.utils import Astro_Coords + +class BField(): + def __init__(self,year=None,doy=None,site=1,heights=None,alpha_i=90): + """ + BField class creates an object to get the Magnetic field for a specific date and + height(s). + + Parameters + ---------- + year = A scalar giving the desired year. If the value is None (default value) then + the current year will be used. + doy = A scalar giving the desired day of the year. If the value is None (default va- + lue) then the current doy will be used. + site = An integer to choose the geographic coordinates of the place where the magne- + tic field will be computed. The default value is over Jicamarca (site=1) + heights = An array giving the heights (km) where the magnetic field will be modeled By default the magnetic field will be computed at 100, 500 and 1000km. + alpha_i = Angle to interpolate the magnetic field. + + Modification History + -------------------- + Converted to Object-oriented Programming by Freddy Galindo, ROJ, 07 October 2009. + Added to signal Chain by Joab Apaza, ROJ, Jun 2023. + """ + + tmp = time.localtime() + if year==None: year = tmp[0] + if doy==None: doy = tmp[7] + self.year = year + self.doy = doy + self.site = site + if heights is None: + heights = numpy.array([100,500,1000]) + else: + heights = numpy.array(heights) + self.heights = heights + self.alpha_i = alpha_i + + def getBField(self,maglimits=numpy.array([-37,-37,37,37])): + """ + getBField models the magnetic field for a different heights in a specific date. + + Parameters + ---------- + maglimits = An 4-elements array giving ..... The default value is [-7,-7,7,7]. + + Return + ------ + dcos = An 4-dimensional array giving the directional cosines of the magnetic field + over the desired place. + alpha = An 3-dimensional array giving the angle of the magnetic field over the desi- + red place. + + Modification History + -------------------- + Converted to Python by Freddy R. Galindo, ROJ, 07 October 2009. + """ + + x_ant = numpy.array([1,0,0]) + y_ant = numpy.array([0,1,0]) + z_ant = numpy.array([0,0,1]) + + if self.site==0: + title_site = "Magnetic equator" + coord_site = numpy.array([-76+52./60.,-11+57/60.,0.5]) + elif self.site==1: + title_site = 'Jicamarca' + coord_site = [-76-52./60.,-11-57/60.,0.5] + heta = (45+5.35)*numpy.pi/180. # (50.35 and 1.46 from Fleish Thesis) + delta = -1.46*numpy.pi/180 + + + x_ant1 = numpy.roll(self.rotvector(self.rotvector(x_ant,1,delta),3,theta),1) + y_ant1 = numpy.roll(self.rotvector(self.rotvector(y_ant,1,delta),3,theta),1) + z_ant1 = numpy.roll(self.rotvector(self.rotvector(z_ant,1,delta),3,theta),1) + + ang0 = -1*coord_site[0]*numpy.pi/180. + ang1 = coord_site[1]*numpy.pi/180. + x_ant = self.rotvector(self.rotvector(x_ant1,2,ang1),3,ang0) + y_ant = self.rotvector(self.rotvector(y_ant1,2,ang1),3,ang0) + z_ant = self.rotvector(self.rotvector(z_ant1,2,ang1),3,ang0) + + elif self.site==2: #AMISR + title_site = 'AMISR 14' + + coord_site = [-76.874913, -11.953371, 0.52984] + + theta = (0.0977)*numpy.pi/180. # 0.0977 + delta = 0.110*numpy.pi/180 # 0.11 + + x_ant1 = numpy.roll(self.rotvector(self.rotvector(x_ant,1,delta),3,theta),1) + y_ant1 = numpy.roll(self.rotvector(self.rotvector(y_ant,1,delta),3,theta),1) + z_ant1 = numpy.roll(self.rotvector(self.rotvector(z_ant,1,delta),3,theta),1) + + ang0 = -1*coord_site[0]*numpy.pi/180. + ang1 = coord_site[1]*numpy.pi/180. + x_ant = self.rotvector(self.rotvector(x_ant1,2,ang1),3,ang0) + y_ant = self.rotvector(self.rotvector(y_ant1,2,ang1),3,ang0) + z_ant = self.rotvector(self.rotvector(z_ant1,2,ang1),3,ang0) + else: +# print "No defined Site. Skip..." + return None + + nhei = self.heights.size + pt_intercep = numpy.zeros((nhei,2)) + nfields = 1 + + grid_res = 2.5 + nlon = int(int(maglimits[2] - maglimits[0])/grid_res + 1) + nlat = int(int(maglimits[3] - maglimits[1])/grid_res + 1) + + location = numpy.zeros((nlon,nlat,2)) + mlon = numpy.atleast_2d(numpy.arange(nlon)*grid_res + maglimits[0]) + mrep = numpy.atleast_2d(numpy.zeros(nlat) + 1) + location0 = numpy.dot(mlon.transpose(),mrep) + + mlat = numpy.atleast_2d(numpy.arange(nlat)*grid_res + maglimits[1]) + mrep = numpy.atleast_2d(numpy.zeros(nlon) + 1) + location1 = numpy.dot(mrep.transpose(),mlat) + + location[:,:,0] = location0 + location[:,:,1] = location1 + + alpha = numpy.zeros((nlon,nlat,nhei)) + rr = numpy.zeros((nlon,nlat,nhei,3)) + dcos = numpy.zeros((nlon,nlat,nhei,2)) + + global first_time + + first_time = None + for ilon in numpy.arange(nlon): + for ilat in numpy.arange(nlat): + outs = self.__bdotk(self.heights, + self.year + self.doy/366., + coord_site[1], + coord_site[0], + coord_site[2], + coord_site[1]+location[ilon,ilat,1], + location[ilon,ilat,0]*720./180.) + + alpha[ilon, ilat,:] = outs[1] + rr[ilon, ilat,:,:] = outs[3] + + mrep = numpy.atleast_2d((numpy.zeros(nhei)+1)).transpose() + tmp = outs[3]*numpy.dot(mrep,numpy.atleast_2d(x_ant)) + tmp = tmp.sum(axis=1) + dcos[ilon,ilat,:,0] = tmp/numpy.sqrt((outs[3]**2).sum(axis=1)) + + mrep = numpy.atleast_2d((numpy.zeros(nhei)+1)).transpose() + tmp = outs[3]*numpy.dot(mrep,numpy.atleast_2d(y_ant)) + tmp = tmp.sum(axis=1) + dcos[ilon,ilat,:,1] = tmp/numpy.sqrt((outs[3]**2).sum(axis=1)) + + return dcos, alpha, nlon, nlat + + + def __bdotk(self,heights,tm,gdlat=-11.95,gdlon=-76.8667,gdalt=0.0,decd=-12.88, ham=-4.61666667): + + global first_time + # Mean Earth radius in Km WGS 84 + a_igrf = 6371.2 + + bk = numpy.zeros(heights.size) + alpha = numpy.zeros(heights.size) + bfm = numpy.zeros(heights.size) + rr = numpy.zeros((heights.size,3)) + rgc = numpy.zeros((heights.size,3)) + + ObjGeodetic = Astro_Coords.Geodetic(gdlat,gdalt) + [gclat,gcalt] = ObjGeodetic.change2geocentric() + + gclat = gclat*numpy.pi/180. + gclon = gdlon*numpy.pi/180. + + # Antenna position from center of Earth + ca_vector = [numpy.cos(gclat)*numpy.cos(gclon),numpy.cos(gclat)*numpy.sin(gclon),numpy.sin(gclat)] + ca_vector = gcalt*numpy.array(ca_vector) + + dec = decd*numpy.pi/180. + + # K vector respect to the center of earth. + klon = gclon + ham*numpy.pi/720. + k_vector = [numpy.cos(dec)*numpy.cos(klon),numpy.cos(dec)*numpy.sin(klon),numpy.sin(dec)] + k_vector = numpy.array(k_vector) + + for ih in numpy.arange(heights.size): + # Vector from Earth's center to volume of interest + rr[ih,:] = k_vector*heights[ih] + cv_vector = numpy.squeeze(ca_vector) + rr[ih,:] + + cv_gcalt = numpy.sqrt(numpy.sum(cv_vector**2.)) + cvxy = numpy.sqrt(numpy.sum(cv_vector[0:2]**2.)) + + radial = cv_vector/cv_gcalt + east = numpy.array([-1*cv_vector[1],cv_vector[0],0])/cvxy + comp1 = east[1]*radial[2] - radial[1]*east[2] + comp2 = east[2]*radial[0] - radial[2]*east[0] + comp3 = east[0]*radial[1] - radial[0]*east[1] + north = -1*numpy.array([comp1, comp2, comp3]) + + rr_k = cv_vector - numpy.squeeze(ca_vector) + u_rr = rr_k/numpy.sqrt(numpy.sum(rr_k**2.)) + + cv_gclat = numpy.arctan2(cv_vector[2],cvxy) + cv_gclon = numpy.arctan2(cv_vector[1],cv_vector[0]) + + bhei = cv_gcalt-a_igrf + blat = cv_gclat*180./numpy.pi + blon = cv_gclon*180./numpy.pi + bfield = self.__igrfkudeki(bhei,tm,blat,blon) + + B = (bfield[0]*north + bfield[1]*east - bfield[2]*radial)*1.0e-5 + + bfm[ih] = numpy.sqrt(numpy.sum(B**2.)) #module + bk[ih] = numpy.sum(u_rr*B) + alpha[ih] = numpy.arccos(bk[ih]/bfm[ih])*180/numpy.pi + rgc[ih,:] = numpy.array([cv_gclon, cv_gclat, cv_gcalt]) + + return bk, alpha, bfm, rr, rgc + + + def __igrfkudeki(self,heights,time,latitude,longitude,ae=6371.2): + """ + __igrfkudeki calculates the International Geomagnetic Reference Field for given in- + put conditions based on IGRF2005 coefficients. + + Parameters + ---------- + heights = Scalar or vector giving the height above the Earth of the point in ques- + tion in kilometers. + time = Scalar or vector giving the decimal year of time in question (e.g. 1991.2). + latitude = Latitude of point in question in decimal degrees. Scalar or vector. + longitude = Longitude of point in question in decimal degrees. Scalar or vector. + ae = + first_time = + + Return + ------ + bn = + be = + bd = + bmod = + balpha = + first_time = + + Modification History + -------------------- + Converted to Python by Freddy R. Galindo, ROJ, 03 October 2009. + """ + + global first_time + global gs, hs, nvec, mvec, maxcoef + + heights = numpy.atleast_1d(heights) + time = numpy.atleast_1d(time) + latitude = numpy.atleast_1d(latitude) + longitude = numpy.atleast_1d(longitude) + + if numpy.max(latitude)==90: +# print "Field calculations are not supported at geographic poles" + pass + + # output arrays + bn = numpy.zeros(heights.size) + be = numpy.zeros(heights.size) + bd = numpy.zeros(heights.size) + + if first_time==None:first_time=0 + + time0 = time[0] + if time!=first_time: + #print "Getting coefficients for", time0 + [periods,g,h ] = self.__readIGRFcoeff() + top_year = numpy.max(periods) + nperiod = (top_year - 1900)/5 + 1 + + maxcoef = 10 + + if time0>=2000:maxcoef = 12 + + + # Normalization array for Schmidt fucntions + multer = numpy.zeros((2+maxcoef,1+maxcoef)) + 1 + for cn in (numpy.arange(maxcoef)+1): + for rm in (numpy.arange(cn)+1): + tmp = numpy.arange(2*rm) + cn - rm + 1. + multer[rm+1,cn] = ((-1.)**rm)*numpy.sqrt(2./tmp.prod()) + + schmidt = multer[1:,1:].transpose() + + # n and m arrays + nvec = numpy.atleast_2d(numpy.arange(maxcoef)+2) + mvec = numpy.atleast_2d(numpy.arange(maxcoef+1)).transpose() + + # Time adjusted igrf g and h with Schmidt normalization + # IGRF coefficient arrays: g0(n,m), n=1, maxcoeff,m=0, maxcoeff, ... + if time0 top_year + dtime = (time0 - top_year) + 5 + ntime = int(g[:,0,0].size - 2) + + + g0 = g[ntime,1:maxcoef+1,:maxcoef+1] + h0 = h[ntime,1:maxcoef+1,:maxcoef+1] + gdot = g[ntime+1,1:maxcoef+1,:maxcoef+1]-g[ntime,1:maxcoef+1,:maxcoef+1] + hdot = h[ntime+1,1:maxcoef+1,:maxcoef+1]-h[ntime,1:maxcoef+1,:maxcoef+1] + gs = (g0 + dtime*(gdot/5.))*schmidt[:maxcoef,0:maxcoef+1] + hs = (h0 + dtime*(hdot/5.))*schmidt[:maxcoef,0:maxcoef+1] + + first_time = time0 + + for ii in numpy.arange(heights.size): + # Height dependence array rad = (ae/(ae+height))**(n+3) + rad = numpy.atleast_2d((ae/(ae + heights[ii]))**(nvec+1)) + + # Sin and Cos of m times longitude phi arrays + mphi = mvec*longitude[ii]*numpy.pi/180. + cosmphi = numpy.atleast_2d(numpy.cos(mphi)) + sinmphi = numpy.atleast_2d(numpy.sin(mphi)) + + # Cos of colatitude theta + c = numpy.cos((90 - latitude[ii])*numpy.pi/180.) + + # Legendre functions p(n,m|c) + [p,dp]= lpmn(maxcoef+1,maxcoef+1,c) + p = p[:,:-1].transpose() + s = numpy.sqrt((1. - c)*(1 + c)) + + # Generate derivative array dpdtheta = -s*dpdc + dpdtheta = c*p/s + for m in numpy.arange(maxcoef+2): dpdtheta[:,m] = m*dpdtheta[:,m] + dpdtheta = dpdtheta + numpy.roll(p,-1,axis=1) + + # Extracting arrays required for field calculations + p = p[1:maxcoef+1,:maxcoef+1] + dpdtheta = dpdtheta[1:maxcoef+1,:maxcoef+1] + + # Weigh p and dpdtheta with gs and hs coefficients. + gp = gs*p + hp = hs*p + gdpdtheta = gs*dpdtheta + hdpdtheta = hs*dpdtheta + # Calcultate field components + matrix0 = numpy.dot(gdpdtheta,cosmphi) + matrix1 = numpy.dot(hdpdtheta,sinmphi) + bn[ii] = numpy.dot(rad,(matrix0 + matrix1)) + matrix0 = numpy.dot(hp,(mvec*cosmphi)) + matrix1 = numpy.dot(gp,(mvec*sinmphi)) + be[ii] = numpy.dot((-1*rad),((matrix0 - matrix1)/s)) + matrix0 = numpy.dot(gp,cosmphi) + matrix1 = numpy.dot(hp,sinmphi) + bd[ii] = numpy.dot((-1*nvec*rad),(matrix0 + matrix1)) + + bmod = numpy.sqrt(bn**2. + be**2. + bd**2.) + btheta = numpy.arctan(bd/numpy.sqrt(be**2. + bn**2.))*180/numpy.pi + balpha = numpy.arctan(be/bn)*180./numpy.pi + + #bn : north + #be : east + #bn : radial + #bmod : module + + + return bn, be, bd, bmod, btheta, balpha + + def str2num(self, datum): + try: + return int(datum) + except: + try: + return float(datum) + except: + return datum + + def __readIGRFfile(self, filename): + list_years=[] + for i in range(1,26): + list_years.append(1895.0 + i*5) + + epochs=list_years + epochs.append(epochs[-1]+5) + nepochs = numpy.shape(epochs) + + gg = numpy.zeros((13,14,nepochs[0]),dtype=float) + hh = numpy.zeros((13,14,nepochs[0]),dtype=float) + + coeffs_file=open(filename) + lines=coeffs_file.readlines() + + coeffs_file.close() + + for line in lines: + items = line.split() + g_h = items[0] + n = self.str2num(items[1]) + m = self.str2num(items[2]) + + coeffs = items[3:] + + for i in range(len(coeffs)-1): + coeffs[i] = self.str2num(coeffs[i]) + + #coeffs = numpy.array(coeffs) + ncoeffs = numpy.shape(coeffs)[0] + + if g_h == 'g': + # print n," g ",m + gg[n-1,m,:]=coeffs + elif g_h=='h': + # print n," h ",m + hh[n-1,m,:]=coeffs + # else : + # continue + + # Ultimo Reordenamiento para almacenar . + gg[:,:,nepochs[0]-1] = gg[:,:,nepochs[0]-2] + 5*gg[:,:,nepochs[0]-1] + hh[:,:,nepochs[0]-1] = hh[:,:,nepochs[0]-2] + 5*hh[:,:,nepochs[0]-1] + +# return numpy.array([gg,hh]) + periods = numpy.array(epochs) + g = gg + h = hh + return periods, g, h + + + def __readIGRFcoeff(self,filename="igrf10coeffs.dat"): + """ + __readIGRFcoeff reads the coefficients from a binary file which is located in the + folder "resource." + + Parameter + --------- + filename = A string to specify the name of the file which contains thec coeffs. The + default value is "igrf10coeffs.dat" + + Return + ------ + periods = A lineal array giving... + g1 = + h1 = + + Modification History + -------------------- + Converted to Python by Freddy R. Galindo, ROJ, 03 October 2009. + """ + + base_path = os.path.dirname(os.path.abspath(__file__)) + filename = os.path.join(base_path, "igrf13coeffs.txt") + + period_v, g_v, h_v = self.__readIGRFfile(filename) + g2 = numpy.zeros((14,14,26)) + h2 = numpy.zeros((14,14,26)) + g2[1:14,:,:] = g_v + h2[1:14,:,:] = h_v + + g = numpy.transpose(g2, (2,0,1)) + h = numpy.transpose(h2, (2,0,1)) + periods = period_v.copy() + + return periods, g, h + + def rotvector(self,vector,axis=1,ang=0): + """ + rotvector function returns the new vector generated rotating the rectagular coords. + + Parameters + ---------- + vector = A lineal 3-elements array (x,y,z). + axis = A integer to specify the axis used to rotate the coord systems. The default + value is 1. + axis = 1 -> Around "x" + axis = 2 -> Around "y" + axis = 3 -> Around "z" + ang = Angle of rotation (in radians). The default value is zero. + + Return + ------ + rotvector = A lineal array of 3 elements giving the new coordinates. + + Modification History + -------------------- + Converted to Python by Freddy R. Galindo, ROJ, 01 October 2009. + """ + + if axis==1: + t = [[1,0,0],[0,numpy.cos(ang),numpy.sin(ang)],[0,-numpy.sin(ang),numpy.cos(ang)]] + elif axis==2: + t = [[numpy.cos(ang),0,-numpy.sin(ang)],[0,1,0],[numpy.sin(ang),0,numpy.cos(ang)]] + elif axis==3: + t = [[numpy.cos(ang),numpy.sin(ang),0],[-numpy.sin(ang),numpy.cos(ang),0],[0,0,1]] + + rotvector = numpy.array(numpy.dot(numpy.array(t),numpy.array(vector))) + + return rotvector diff --git a/schainpy/model/utils/Misc_Routines.py b/schainpy/model/utils/Misc_Routines.py new file mode 100644 index 0000000..bc398e3 --- /dev/null +++ b/schainpy/model/utils/Misc_Routines.py @@ -0,0 +1,61 @@ +""" +The module MISC_ROUTINES gathers classes and functions which are useful for daily processing. As an +example we have conversion factor or universal constants. + +MODULES CALLED: +NUMPY, SYS + +MODIFICATION HISTORY: +Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ, 21 October 2009. +""" + +import numpy +import sys + +class CoFactors(): + """ + CoFactor class used to call pre-defined conversion factor (e.g. degree to radian). The cu- + The current available factor are: + + d2r = degree to radian. + s2r = seconds to radian?, degree to arcsecond.? + h2r = hour to radian. + h2d = hour to degree + """ + + d2r = numpy.pi/180. + s2r = numpy.pi/(180.*3600.) + h2r = numpy.pi/12. + h2d = 15. + + +class Vector: + """ + direction = 0 Polar to rectangular; direction=1 rectangular to polar + """ + def __init__(self,vect,direction=0): + nsize = numpy.size(vect) + if nsize <= 3: + vect = vect.reshape(1,nsize) + + self.vect = vect + self.dirc = direction + + + + def Polar2Rect(self): + if self.dirc == 0: + jvect = self.vect*numpy.pi/180. + mmx = numpy.cos(jvect[:,1])*numpy.sin(jvect[:,0]) + mmy = numpy.cos(jvect[:,1])*numpy.cos(jvect[:,0]) + mmz = numpy.sin(jvect[:,1]) + mm = numpy.array([mmx,mmy,mmz]).transpose() + + elif self.dirc == 1: + mm = [numpy.arctan2(self.vect[:,0],self.vect[:,1]),numpy.arcsin(self.vect[:,2])] + mm = numpy.array(mm)*180./numpy.pi + + return mm + + + \ No newline at end of file diff --git a/schainpy/model/utils/TimeTools.py b/schainpy/model/utils/TimeTools.py new file mode 100644 index 0000000..1750f48 --- /dev/null +++ b/schainpy/model/utils/TimeTools.py @@ -0,0 +1,430 @@ +""" +The TIME_CONVERSIONS.py module gathers classes and functions for time system transformations +(e.g. between seconds from 1970 to datetime format). + +MODULES CALLED: +NUMPY, TIME, DATETIME, CALENDAR + +MODIFICATION HISTORY: +Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ Aug 13, 2009. +""" + +import numpy +import time +from datetime import datetime as dt +import calendar + +class Time: + """ + time(year,month,dom,hour,min,secs) + + An object represents a date and time of certain event.. + + Parameters + ---------- + YEAR = Number of the desired year. Year must be valid values from the civil calendar. + Years B.C.E must be represented as negative integers. Years in the common era are repre- + sented as positive integers. In particular, note that there is no year 0 in the civil + calendar. 1 B.C.E. (-1) is followed by 1 C.E. (1). + + MONTH = Number of desired month (1=Jan, ..., 12=December). + + DOM = Number of day of the month. + + HOUR = Number of the hour of the day. By default hour=0 + + MINS = Number of the minute of the hour. By default min=0 + + SECS = Number of the second of the minute. By default secs=0. + + Examples + -------- + time_info = time(2008,9,30,12,30,00) + + time_info = time(2008,9,30) + """ + + def __init__(self, year=None, month=None, dom=None, hour=0, mins=0, secs=0): + # If one the first three inputs are not defined, it takes the current date. + date = time.localtime() + if year==None:year=date[0] + if month==None:month=date[1] + if dom==None:dom=date[2] + + # Converting to arrays + year = numpy.array([year]); month = numpy.array([month]); dom = numpy.array([dom]) + hour = numpy.array([hour]); mins = numpy.array([mins]); secs = numpy.array([secs]) + + # Defining time information object. + self.year = numpy.atleast_1d(year) + self.month = numpy.atleast_1d(month) + self.dom = numpy.atleast_1d(dom) + self.hour = numpy.atleast_1d(hour) + self.mins = numpy.atleast_1d(mins) + self.secs = numpy.atleast_1d(secs) + + def change2julday(self): + """ + Converts a datetime to Julian days. + """ + + # Defining constants + greg = 2299171 # incorrect Julian day for Oct, 25, 1582. + min_calendar = -4716 + max_calendar = 5000000 + + min_year = numpy.nanmin(self.year) + max_year = numpy.nanmax(self.year) + if (min_yearmax_calendar): + print ("Value of Julian date is out of allowed range") + return -1 + + noyear = numpy.sum(self.year==0) + if noyear>0: + print ("There is no year zero in the civil calendar") + return -1 + + # Knowing if the year is less than 0. + bc = self.year<0 + + # Knowing if the month is less than March. + inJanFeb = self.month<=2 + + jy = self.year + bc - inJanFeb + jm = self.month + (1 + 12*inJanFeb) + + # Computing Julian days. + jul= numpy.floor(365.25*jy) + numpy.floor(30.6001*jm) + (self.dom+1720995.0) + + # Test whether to change to Gregorian Calendar + if numpy.min(jul) >= greg: + ja = numpy.int32(0.01*jy) + jul = jul + 2 - ja + numpy.int32(0.25*ja) + else: + gregchange = numpy.where(jul >= greg) + if gregchange[0].size>0: + ja = numpy.int32(0.01 + jy[gregchange]) + jy[gregchange] = jy[gregchange] + 2 - ja + numpy.int32(0.25*ja) + + # Determining machine-specific parameters affecting floating-point. + eps = 0.0 # Replace this line for a function to get precision. + eps = abs(jul)*0.0 > eps + + jul = jul + (self.hour/24. -0.5) + (self.mins/1440.) + (self.secs/86400.) + eps + + return jul[0] + + def change2secs(self): + """ + Converts datetime to number of seconds respect to 1970. + """ + + dtime = dt(self.year[0], self.month[0], self.dom[0]) + return (dtime-dt(1970, 1, 1)).total_seconds() + + year = self.year + if year.size>1: year = year[0] + + month = self.month + if month.size>1: month = month[0] + + dom = self.dom + if dom.size>1: dom = dom[0] + + # Resizing hour, mins and secs if it was necessary. + hour = self.hour + if hour.size>1:hour = hour[0] + if hour.size==1:hour = numpy.resize(hour,year.size) + + mins = self.mins + if mins.size>1:mins = mins[0] + if mins.size==1:mins = numpy.resize(mins,year.size) + + secs = self.secs + if secs.size>1:secs = secs[0] + if secs.size==1:secs = numpy.resize(secs,year.size) + + # Using time.mktime to compute seconds respect to 1970. + secs1970 = numpy.zeros(year.size) + for ii in numpy.arange(year.size): + secs1970[ii] = time.mktime((int(year[ii]),int(month[ii]),int(dom[ii]),\ + int(hour[ii]),int(mins[ii]),int(secs[ii]),0,0,0)) + + secs1970 = numpy.int32(secs1970 - time.timezone) + + return secs1970 + + def change2strdate(self,mode=1): + """ + change2strdate method converts a date and time of certain event to date string. The + string format is like localtime (e.g. Fri Oct 9 15:00:19 2009). + + Parameters + ---------- + None. + + Return + ------ + + Modification History + -------------------- + Created by Freddy R. Galindo, ROJ, 09 October 2009. + + """ + + secs = numpy.atleast_1d(self.change2secs()) + strdate = [] + for ii in numpy.arange(numpy.size(secs)): + secs_tmp = time.localtime(secs[ii] + time.timezone) + if mode==1: + strdate.append(time.strftime("%d-%b-%Y (%j) %H:%M:%S",secs_tmp)) + elif mode==2: + strdate.append(time.strftime("%d-%b-%Y (%j)",secs_tmp)) + + strdate = numpy.array(strdate) + + return strdate + + +class Secs: + """ + secs(secs): + + An object represents the number of seconds respect to 1970. + + Parameters + ---------- + + SECS = A scalar or array giving the number of seconds respect to 1970. + + Example: + -------- + secs_info = secs(1251241373) + + secs_info = secs([1251241373,1251241383,1251241393]) + """ + def __init__(self,secs): + self.secs = secs + + def change2julday(self): + """ + Convert seconds from 1970 to Julian days. + """ + + secs_1970 = time(1970,1,1,0,0,0).change2julday() + + julian = self.secs/86400.0 + secs_1970 + + return julian + + def change2time(self): + """ + Converts seconds from 1970 to datetime. + """ + + secs1970 = numpy.atleast_1d(self.secs) + + datetime = numpy.zeros((9,secs1970.size)) + for ii in numpy.arange(secs1970.size): + tuple = time.gmtime(secs1970[ii]) + datetime[0,ii] = tuple[0] + datetime[1,ii] = tuple[1] + datetime[2,ii] = tuple[2] + datetime[3,ii] = tuple[3] + datetime[4,ii] = tuple[4] + datetime[5,ii] = tuple[5] + datetime[6,ii] = tuple[6] + datetime[7,ii] = tuple[7] + datetime[8,ii] = tuple[8] + + datetime = numpy.int32(datetime) + + return datetime + + +class Julian: + """ + julian(julian): + + An object represents julian days. + + Parameters + ---------- + + JULIAN = A scalar or array giving the julina days. + + Example: + -------- + julian_info = julian(2454740) + + julian_info = julian([2454740,2454760,2454780]) + """ + def __init__(self,julian): + self.julian = numpy.atleast_1d(julian) + + def change2time(self): + """ + change2time method converts from julian day to calendar date and time. + + Return + ------ + year = An array giving the year of the desired julian day. + month = An array giving the month of the desired julian day. + dom = An array giving the day of the desired julian day. + hour = An array giving the hour of the desired julian day. + mins = An array giving the minute of the desired julian day. + secs = An array giving the second of the desired julian day. + + Examples + -------- + >> jd = 2455119.0 + >> [yy,mo,dd,hh,mi,ss] = TimeTools.julian(jd).change2time() + >> print [yy,mo,dd,hh,mi,ss] + [2009] [10] [ 14.] [ 12.] [ 0.] [ 0.] + + Modification history + -------------------- + Translated from "Numerical Recipies in C", by William H. Press, Brian P. Flannery, + Saul A. Teukolsky, and William T. Vetterling. Cambridge University Press, 1988. + Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. + """ + + min_julian = -1095 + max_julian = 1827933925 + if (numpy.min(self.julian) < min_julian) or (numpy.max(self.julian) > max_julian): + print ('Value of Julian date is out of allowed range.') + return None + + # Beginning of Gregorian calendar + igreg = 2299161 + julLong = numpy.floor(self.julian + 0.5) + minJul = numpy.min(julLong) + + if (minJul >= igreg): + # All are Gregorian + jalpha = numpy.int32(((julLong - 1867216) - 0.25)/36524.25) + ja = julLong + 1 + jalpha - numpy.int32(0.25*jalpha) + else: + ja = julLong + gregChange = numpy.where(julLong >= igreg) + if gregChange[0].size>0: + jalpha = numpy.int32(((julLong[gregChange]-1867216) - 0.25)/36524.25) + ja[gregChange] = julLong[gregChange]+1+jalpha-numpy.int32(0.25*jalpha) + + # clear memory. + jalpha = -1 + + jb = ja + 1524 + jc = numpy.int32(6680. + ((jb-2439870)-122.1)/365.25) + jd = numpy.int32(365.*jc + (0.25*jc)) + je = numpy.int32((jb - jd)/30.6001) + + dom = jb - jd - numpy.int32(30.6001*je) + month = je - 1 + month = ((month - 1) % 12) + 1 + month = numpy.atleast_1d(month) + year = jc - 4715 + year = year - (month > 2)*1 + year = year - (year <= 0)*1 + year = numpy.atleast_1d(year) + + # Getting hours, minutes, seconds + fraction = self.julian + 0.5 - julLong + eps_0 = dom*0.0 + 1.0e-12 + eps_1 = 1.0e-12*numpy.abs(julLong) + eps = (eps_0>eps_1)*eps_0 + (eps_0<=eps_1)*eps_1 + + hour_0 = dom*0 + 23 + hour_2 = dom*0 + 0 + hour_1 = numpy.floor(fraction*24.0 + eps) + hour = ((hour_1>hour_0)*23) + ((hour_1<=hour_0)*hour_1) + hour = ((hour_1=hour_2)*hour_1) + + fraction = fraction - (hour/24.0) + mins_0 = dom*0 + 59 + mins_2 = dom*0 + 0 + mins_1 = numpy.floor(fraction*1440.0 + eps) + mins = ((mins_1>mins_0)*59) + ((mins_1<=mins_0)*mins_1) + mins = ((mins_1=mins_2)*mins_1) + + secs_2 = dom*0 + 0 + secs_1 = (fraction - mins/1440.0)*86400.0 + secs = ((secs_1=secs_2)*secs_1) + + return year,month,dom,hour,mins,secs + + def change2secs(self): + """ + Converts from Julian days to seconds from 1970. + """ + + jul_1970 = Time(1970,1,1,0,0,0).change2julday() + + secs = numpy.int32((self.julian - jul_1970)*86400) + + return secs + + def change2lst(self,longitude=-76.874369): + """ + CT2LST converts from local civil time to local mean sideral time + + longitude = The longitude in degrees (east of Greenwich) of the place for which + the local sideral time is desired, scalar. The Greenwich mean sideral time (GMST) + can be found by setting longitude=0. + """ + + # Useful constants, see Meus, p. 84 + c = numpy.array([280.46061837, 360.98564736629, 0.000387933, 38710000.0]) + jd2000 = 2451545.0 + t0 = self.julian - jd2000 + t = t0/36525. + + # Computing GST in seconds + theta = c[0] + (c[1]*t0) + (t**2)*(c[2]-t/c[3]) + + # Computing LST in hours + lst = (theta + longitude)/15.0 + neg = numpy.where(lst < 0.0) + if neg[0].size>0:lst[neg] = 24.0 + (lst[neg] % 24) + lst = lst % 24.0 + + return lst + + +class date2doy: + def __init__(self,year,month,day): + self.year = year + self.month = month + self.day = day + + def change2doy(self): + if calendar.isleap(self.year) == True: + tfactor = 1 + else: + tfactor = 2 + + day = self.day + month = self.month + + doy = numpy.floor((275*month)/9.0) - (tfactor*numpy.floor((month+9)/12.0)) + day - 30 + + return numpy.int32(doy) + + +class Doy2Date: + def __init__(self,year,doy): + self.year = year + self.doy = doy + + def change2date(self): + months = numpy.arange(12) + 1 + + first_dem = date2doy(self.year,months,1) + first_dem = first_dem.change2doy() + + imm = numpy.where((self.doy - first_dem) > 0) + + month = imm[0].size + dom = self.doy -first_dem[month - 1] + 1 + + return month, dom