1 """
2 Calculate stresses on an icy satellite over a rectangular geographic region on
3 a regularly spaced lat-lon grid.
4
5 The datacube containing the results of the calculation are output as a
6 U{Unidata netCDF <http://www.unidata.ucar.edu/software/netcdf>} (.nc) file,
7 which can be displayed using a wide variety of visualization software.
8
9 """
10
11 import SatStress as SS
12 import time
13 import netCDF3
14 import physcon as pc
15 import numpy
16 from numpy import linalg
17 import scipy
18 from optparse import OptionParser
19
20
22 """
23 Calculate satellite stresses on a regular grid within a lat-lon window.
24
25 """
26
27 usage = "usage: %prog [options] satfile gridfile outfile"
28 description = __doc__
29 op = OptionParser(usage)
30
31 (options, args) = op.parse_args()
32
33
34 if len(args) != 3:
35 op.error("incorrect number of arguments")
36
37
38 satfile = open(args[0], 'r')
39 the_sat = SS.Satellite(satfile)
40 satfile.close()
41
42
43 gridfile = open(args[1], 'r')
44 the_grid = Grid(gridfile, satellite=the_sat)
45 gridfile.close()
46
47
48 diurnal_stress = SS.StressCalc([SS.Diurnal(the_sat),])
49 nsr_stress = SS.StressCalc([SS.NSR(the_sat),])
50
51
52
53 nc_out = netCDF3.Dataset(args[2], 'w')
54
55
56
57
58
59
60
61
62
63
64
65
66 nc_out.description = "Testing pySatStress on a regular grid"
67 nc_out.history = "Created: %s using pySatStress" % ( time.ctime(time.time()) )
68 nc_out.Conventions = "COARDS"
69
70
71
72 n_deltas = the_grid.nsr_delta_numsteps
73 nc_out.createDimension('nsr_delta', n_deltas)
74 nsr_deltas = nc_out.createVariable('nsr_delta', 'f4', ('nsr_delta',))
75 nsr_deltas.units = ""
76 nsr_deltas.long_name = "NSR Surface Delta (mu/(eta*omega))"
77 nsr_deltas[:] = numpy.logspace(the_grid.nsr_delta_min, the_grid.nsr_delta_max, num=n_deltas)
78
79 n_times = len(numpy.arange(the_grid.time_min, the_grid.time_max, the_grid.time_step))
80 nc_out.createDimension('time', n_times)
81 times = nc_out.createVariable('time', 'f4', ('time',))
82
83
84
85 if the_grid.orbit_min is None:
86 times.units = "seconds"
87 times.long_name = "time after periapse"
88 times[:] = numpy.arange(the_grid.time_min, the_grid.time_max, the_grid.time_step)
89 else:
90 times.units = "degrees"
91 times.long_name = "degrees after periapse"
92 times[:] = numpy.arange(the_grid.orbit_min, the_grid.orbit_max, the_grid.orbit_step)
93
94 n_lats = len(numpy.arange(the_grid.lat_min, the_grid.lat_max+the_grid.latlon_step, the_grid.latlon_step))
95 nc_out.createDimension('latitude', n_lats)
96 lats = nc_out.createVariable('latitude', 'f4', ('latitude',))
97 lats.units = "degrees_north"
98 lats.long_name = "latitude"
99 lats[:] = numpy.arange(the_grid.lat_min, the_grid.lat_max + (the_grid.latlon_step/2.0), the_grid.latlon_step)
100
101 n_lons = len(numpy.arange(the_grid.lon_min, the_grid.lon_max+the_grid.latlon_step, the_grid.latlon_step))
102 nc_out.createDimension('longitude', n_lons)
103 lons = nc_out.createVariable('longitude', 'f4', ('longitude',))
104 lons.units = "degrees_east"
105 lons.long_name = "longitude"
106 lons[:] = numpy.arange(the_grid.lon_min, the_grid.lon_max + (the_grid.latlon_step/2.0), the_grid.latlon_step)
107
108
109
110 Ttt_D = nc_out.createVariable('Ttt_D', 'f4', ('time', 'latitude', 'longitude',))
111 Ttt_D.units = "Pa"
112 Ttt_D.long_name = "north-south component stress of Diurnal stresses"
113
114 Tpt_D = nc_out.createVariable('Tpt_D', 'f4', ('time', 'latitude', 'longitude',))
115 Tpt_D.units = "Pa"
116 Tpt_D.long_name = "shear component of Diurnal stresses"
117
118 Tpp_D = nc_out.createVariable('Tpp_D', 'f4', ('time', 'latitude', 'longitude',))
119 Tpp_D.units = "Pa"
120 Tpp_D.long_name = "east-west component of Diurnal stresses"
121
122 Ttt_N = nc_out.createVariable('Ttt_N', 'f4', ('nsr_delta', 'latitude', 'longitude',))
123 Ttt_N.units = "Pa"
124 Ttt_N.long_name = "north-south component of NSR stresses"
125
126 Tpt_N = nc_out.createVariable('Tpt_N', 'f4', ('nsr_delta', 'latitude', 'longitude',))
127 Tpt_N.units = "Pa"
128 Tpt_N.long_name = "shear component of NSR stresses"
129
130 Tpp_N = nc_out.createVariable('Tpp_N', 'f4', ('nsr_delta', 'latitude', 'longitude',))
131 Tpp_N.units = "Pa"
132 Tpp_N.long_name = "east-west component of NSR stresses"
133
134
135
136 for t in range(n_times):
137
138
139
140
141 if the_grid.orbit_min is None:
142 time_sec = times[t]
143 else:
144 time_sec = the_sat.orbit_period()*(times[t]/360.0)
145
146 print "Calculating diurnal stresses at", times[t], times.long_name
147 for lon in range(n_lons):
148 for lat in range(n_lats):
149
150 Tau_D = diurnal_stress.tensor(theta = scipy.radians(90.0-lats[lat]),\
151 phi = scipy.radians(lons[lon]),\
152 t = time_sec )
153
154 nc_out.variables['Ttt_D'][t,lat,lon] = Tau_D[0,0]
155 nc_out.variables['Tpt_D'][t,lat,lon] = Tau_D[1,0]
156 nc_out.variables['Tpp_D'][t,lat,lon] = Tau_D[1,1]
157
158
159
160
161
162 nc_out.sync()
163
164
165
166 for delta in range(n_deltas):
167
168 the_sat.nsr_period = 4*numpy.pi*nsr_deltas[delta]*the_sat.layers[-1].viscosity/the_sat.layers[-1].lame_mu
169
170 nsr_stress = SS.StressCalc([SS.NSR(the_sat),])
171
172 print "Calculating NSR stresses at Delta =", nsr_deltas[delta], "(P_nsr = %g yr)" % (the_sat.nsr_period/31556926)
173 for lon in range(n_lons):
174 for lat in range(n_lats):
175
176
177
178
179
180
181
182 Tau_N = nsr_stress.tensor(theta = scipy.radians(90.0-lats[lat]),\
183 phi = scipy.radians(lons[lon]),\
184 t = 0)
185
186 nc_out.variables['Ttt_N'][delta,lat,lon] = Tau_N[0,0]
187 nc_out.variables['Tpt_N'][delta,lat,lon] = Tau_N[1,0]
188 nc_out.variables['Tpp_N'][delta,lat,lon] = Tau_N[1,1]
189
190
191
192
193
194
195
196
197
198
199 nc_out.sync()
200
201
202
204 """
205 A container class defining the temporal and geographic range and resolution
206 of the calculation.
207
208 The parameters defining the calculation grid are read in from a name value
209 file, parsed into a Python dictionary using L{SatStress.nvf2dict}, and used
210 to set the data attributes of the L{Grid} object.
211
212 The geographic extent of the calculation is specified by minimum and maximum
213 values for latitude and longitude.
214
215 The geographic resolution of the calculation is defined by an angular
216 separation between calculations. This angular separation is the same in
217 the north-south and the east-west direction.
218
219 The temporal range and resolution of the calculation can be specified either
220 in terms of actual time units (seconds) or in terms of the satellite's orbital
221 position (in degrees). In both cases, time=0 is taken to occur at periapse.
222
223 S{Delta} is a measure of how viscous or elastic the response of the body
224 is. It's equal to (S{mu})/(S{eta}S{omega}) where S{mu} and S{eta} are the
225 shear modulus and viscosity of the surface layer, respectively, and
226 S{omega} is the forcing frequency to which the body is subjected (see Wahr
227 et al. (2008) for a detailed discussion). It is a logarithmic parameter, so
228 its bounds are specified as powers of 10, e.g. if the minimum value is -3,
229 the initial S{Delta} is 10^-3 = 0.001.
230
231 @ivar grid_id: A string identifying the grid
232 @type grid_id: str
233 @ivar lat_min: Southern bound, degrees (north positive).
234 @type lat_min: float
235 @ivar lat_max: Northern bound, degrees (north positive).
236 @type lat_max: float
237 @ivar lon_min: Western bound, degrees (east positive).
238 @type lon_min: float
239 @ivar lon_max: Eastern bound, degrees (east positive).
240 @type lon_max: float
241 @ivar latlon_step: Angular separation between calculations.
242 @type latlon_step: float
243 @ivar time_min: Initial time at which calculation begins (0 = periapse).
244 @type time_min: float
245 @ivar time_max: Final time at which calculation ends.
246 @type time_max: float
247 @ivar time_step: Seconds between subsequent calculations.
248 @type time_step: float
249 @ivar orbit_min: Initial orbital position in degrees (0 = periapse)
250 @type orbit_min: float
251 @ivar orbit_max: Final orbital position in degrees (0 = periapse)
252 @type orbit_max: float
253 @ivar orbit_step: Orbital angular separation between calculations in degrees
254 @type orbit_step: float
255 @ivar nsr_delta_min: Initial S{Delta} = 10^(nsr_delta_min)
256 @type nsr_delta_min: float
257 @ivar nsr_delta_max: Final S{Delta} = 10^(nsr_delta_max)
258 @type nsr_delta_max: float
259 @ivar nsr_delta_numsteps: How many S{Delta} values to calculate total
260 @type nsr_delta_numsteps: int
261
262 """
263
264 - def __init__(self, gridFile, satellite=None):
265 """Initialize the Grid object from a gridFile."""
266
267
268 gridParams = SS.nvf2dict(gridFile, comment='#')
269
270 self.grid_id = gridParams['GRID_ID']
271
272 self.lat_min = float(gridParams['LAT_MIN'])
273 self.lat_max = float(gridParams['LAT_MAX'])
274 self.lon_min = float(gridParams['LON_MIN'])
275 self.lon_max = float(gridParams['LON_MAX'])
276 self.latlon_step = float(gridParams['LATLON_STEP'])
277
278 self.nsr_delta_numsteps = float(gridParams['NSR_DELTA_NUMSTEPS'])
279 self.nsr_delta_min = float(gridParams['NSR_DELTA_MIN'])
280 self.nsr_delta_max = float(gridParams['NSR_DELTA_MAX'])
281
282
283 if gridParams.has_key('TIME_MIN'):
284 self.time_min = float(gridParams['TIME_MIN'])
285 self.time_max = float(gridParams['TIME_MAX'])
286 self.time_step = float(gridParams['TIME_STEP'])
287
288 self.orbit_min = None
289 self.orbit_max = None
290 self.orbit_step = None
291
292 elif gridParams.has_key('ORBIT_MIN'):
293 self.orbit_min = float(gridParams['ORBIT_MIN'])
294 self.orbit_max = float(gridParams['ORBIT_MAX'])
295 self.orbit_step = float(gridParams['ORBIT_STEP'])
296
297 self.time_min = satellite.orbit_period()*(self.orbit_min/360.0)
298 self.time_max = satellite.orbit_period()*(self.orbit_max/360.0)
299 self.time_step = satellite.orbit_period()*(self.orbit_step/360.0)
300
301 else:
302
303 pass
304
305
306
307 if __name__ == "__main__":
308 main()
309