Changeset 934 for trunk/UTIL

Timestamp:

Apr 22, 2013, 7:20:23 PM (12 years ago)

Author:

aslmd

Message:

UTIL PYTHON planetoplot_v2. thanks to Thomas Navarro added the possibility to change time axis through attribute self.changetime; for now only mars_sol2ls is available, but easy to add one: put calculations in ppcompute, and add the option in the method performtimechange.

Location:

trunk/UTIL/PYTHON/planetoplot_v2

Files:

• ppclass.py (modified) (6 diffs)
• ppcompute.py (modified) (2 diffs)

trunk/UTIL/PYTHON/planetoplot_v2/ppclass.py

@@ -153 +153 @@
                       color=None,\
                       label=None,\
+                      changetime=None,\
                       title=None):
         self.request = None
@@ -184 +185 @@
         self.folder = folder
         self.includedate = includedate
+        self.changetime = changetime
         ## here are user-defined plot settings 
         ## -- if not None, valid on all plots in the pp() objects
@@ -243 +245 @@
             self.title = other.title
             self.includedate = other.includedate
+            self.changetime = other.changetime
         else:
             print "!! ERROR !! argument must be a pp object." ; exit()
@@ -541 +544 @@
               ### get x,y,z,t dimensions from file
               obj.getdim()
+              ### possible time axis change
+              obj.changetime = self.changetime
+              obj.performtimechange()
               ### get methods 
               obj.method_x = mx ; obj.method_y = my
@@ -1069 +1075 @@
         self.swap_axes = False ; self.invert_axes = False
         self.compute = None
+        self.changetime = None
 
     # open a file. for now it is netcdf. TBD for other formats.
@@ -1181 +1188 @@
           if self.dim_t > 1: 
                if self.verbose: print "**** OK. t axis %4.0f values [%5.1f,%5.1f]" % (self.dim_t,self.field_t.min(),self.field_t.max())     
+
+    # change time axis
+    # ... add your options here!
+    # --------------------------
+    def performtimechange(self):
+        if self.changetime is not None \
+          and self.name_t != "t grid points":
+            if self.verbose: print "**** OK. Converting time axis:",self.changetime
+            if self.changetime == "mars_sol2ls": 
+                self.field_t = ppcompute.mars_sol2ls(self.field_t)
+            else:
+                print "!! WARNING !! This time change is not implemented. Nothing is done."
 
     # get list of index to be retrieved for time axis

trunk/UTIL/PYTHON/planetoplot_v2/ppcompute.py

@@ -7 +7 @@
 # added librairies
 import numpy as np
+import math  as m
 import scipy.signal as sp_signal
 ###############################################
@@ -173 +174 @@
     y=np.convolve(w/w.sum(),s,mode='valid')
     return y
+
+########################
+#### TIME CONVERTER ####
+########################
+
+#  mars_sol2ls
+#  author T. Navarro
+#  -----------------
+#  convert a given martian day number (sol)
+#  into corresponding solar longitude, Ls (in degr.),
+#  where sol=0=Ls=0 is the
+#  northern hemisphere spring equinox.
+#  -----------------
+def mars_sol2ls(soltabin):
+      year_day  = 668.6
+      peri_day  = 485.35
+      e_elips   = 0.09340
+      radtodeg  = 57.2957795130823
+      timeperi  = 1.90258341759902
+      if type(soltabin).__name__ in ['int','float','float32','float64']:
+        soltab=[soltabin]
+        solout=np.zeros([1])
+      else:
+        soltab=soltabin
+        solout=np.zeros([len(soltab)])
+      i=0
+      for sol in soltab:
+         zz=(sol-peri_day)/year_day
+         zanom=2.*np.pi*(zz-np.floor(zz))
+         xref=np.abs(zanom)
+#  The equation zx0 - e * sin (zx0) = xref, solved by Newton
+         zx0=xref+e_elips*m.sin(xref)
+         iter=0
+         while iter <= 10:
+            iter=iter+1
+            zdx=-(zx0-e_elips*m.sin(zx0)-xref)/(1.-e_elips*m.cos(zx0))
+            if(np.abs(zdx) <= (1.e-7)):
+              continue
+            zx0=zx0+zdx
+         zx0=zx0+zdx
+         if(zanom < 0.): zx0=-zx0
+# compute true anomaly zteta, now that eccentric anomaly zx0 is known
+         zteta=2.*m.atan(m.sqrt((1.+e_elips)/(1.-e_elips))*m.tan(zx0/2.))
+# compute Ls
+         ls=zteta-timeperi
+         if(ls < 0.): ls=ls+2.*np.pi
+         if(ls > 2.*np.pi): ls=ls-2.*np.pi
+# convert Ls in deg.
+         ls=radtodeg*ls
+         solout[i]=ls
+         i=i+1
+      return solout
+