Changeset 403 for trunk/UTIL

Timestamp:

Nov 21, 2011, 6:09:07 PM (14 years ago)

Author:

acolaitis

Message:

PYTHON.

############################################
1/------------------------------------------
############################################
# First version of a new command line routine
# that takes in input a diagfi and a MCS datafile
# as formatted by Luca Montabone
# and find the values in the diagfi temp field
# corresponding to daytime and nightime localtimes in MCS file
# (at each longitude and latitude, which is especially important at the poles
# where the spacecraft covers multiple local times for each bin)
# The selected diagfi data is then binned at the Ls of the MCS file
# and written as a diagfi_MCS.nc file.
#
# Example: (you have to be connected on Penn, as the file is on the local disk)

/san0/acolmd/SVN/trunk/UTIL/PYTHON/mcs.py -f /home/local_home/colaitis/Polar_simulations/Polar_All_Th/diagfi16.nc -m /san0/acolmd/DATA/MCS/MCS_processeddata/MCSdata_binned_MY29.nc

# It is then easy to compare the resulting gcm re-arranged data to the MCS dataset using
# your favorite pp command. Note that in this first version of the routine, only the temperature field is
# recomputed. A more complex version will come out later, with variable selection and conservation of
# fields needed to then use zrecast or hrecast.

############################################
2/------------------------------------------
############################################

# Added python versions of ls2sol and sol2ls
# which can work either with single values or arrays
# These routines are in times.py so that you can add
# in your script:

from times import ls2sol,sol2ls

############################################
2/------------------------------------------
############################################

# Added new routine in mymath, which finds the index of the
# element nearest to "value" in "array". The array
# can be masked or not, the output is an array.
# If the value is a NaN, the routine returns a NaN.
# One can specify the option "strict=True" so that
# the routine will return a NaN if the value asked is not
# between the min and max of the array.

Location:

trunk/UTIL/PYTHON

Files:

• mcs.py (added)
• mymath.py (modified) (1 diff)
• time.F (modified) (3 diffs)

trunk/UTIL/PYTHON/mymath.py

@@ -172 +172 @@
     return [item[0] for item in self.items() if item[1] == value]
 
+# A.C. routine that gets the nearest value index of array and value
+def find_nearest(arr,value,axis=None,strict=False):
+    import numpy as np
+    # Special case when the value is nan
+    if value*0 != 0: return np.NaN
+    # Check that the value we search is inside the array for the strict mode
+    if strict:
+       min=arr.min()
+       max=arr.max()
+       if ((value > max) or (value < min)): return np.NaN
+
+    if type(arr).__name__=='MaskedArray':
+       mask=np.ma.getmask(arr)
+       idx=np.ma.argmin(np.abs(arr-value),axis=axis)
+    # Special case when there are only missing values on the axis
+       if mask[idx]:
+          idx=np.NaN
+    else:
+       idx=(np.abs(arr-value)).argmin(axis=axis)
+    return idx
+

trunk/UTIL/PYTHON/time.F

@@ -18 +18 @@
       double precision    pi,radtodeg
 c number of martian days (sols) in a martian year
-      parameter (year_day=668.6d0)
+      parameter (year_day=668.6)
 c perihelion date (in sols)
-      parameter (peri_day=485.35d0)
+      parameter (peri_day=485.35)
 c orbital eccentricity
-      parameter (e_elips=0.09340d0)
-      parameter (pi=3.14159265358979d0)
+      parameter (e_elips=0.09340)
+      parameter (pi=3.14159265358979)
 c  radtodeg: 180/pi
-      parameter (radtodeg=57.2957795130823d0)
+      parameter (radtodeg=57.2957795130823)
 c  timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
-      parameter (timeperi=1.90258341759902d0)
+      parameter (timeperi=1.90258341759902)
 
       double precision    zanom,xref,zx0,zdx,zteta,zz
@@ -41 +41 @@
       do 110 iter=1,10
          zdx=-(zx0-e_elips*dsin(zx0)-xref)/(1.-e_elips*dcos(zx0))
-         if(dabs(zdx).le.(1.d-7)) then ! typically, 2 or 3 iterations are enough
+         if(dabs(zdx).le.(1.d-7)) then
            goto 120
          endif
@@ -77 +77 @@
 c  Local:
       double precision xref,zx0,zteta,zz
-c       xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
-      double precision year_day 
       double precision peri_day,timeperi,e_elips
       double precision pi,degrad 
-      parameter (year_day=668.6d0) ! number of sols in a amartian year
-c      data peri_day /485.0/
-      parameter (peri_day=485.35d0) ! date (in sols) of perihelion
-c  timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
-      parameter (timeperi=1.90258341759902d0)
-      parameter (e_elips=0.0934d0)  ! eccentricity of orbit
-      parameter (pi=3.14159265358979d0)
-      parameter (degrad=57.2957795130823d0)
+      double precision,parameter :: year_day=668.6 
+      parameter (peri_day=485.35) 
+      parameter (timeperi=1.90258341759902)
+      parameter (e_elips=0.0934)  
+      parameter (pi=3.14159265358979)
+      parameter (degrad=57.2957795130823)
 
       if (abs(ls).lt.1.0e-5) then