# Adapted for numpy/ma/cdms2 by convertcdms.py """Documentation for module spheretest: an automatic test for sphere, an interface to spherepack TESTING Typing cdat spheretest.py generates some testing of the spheremodule using analytically generated winds as the input fields. There are 3 tests. test 1 Starting with analytically generated winds, it calculates the stream function and velocity potential using Spherepack and compares with an analytically generated stream function and velocity potential. test 2 Starting with analytically generated winds on a regular grid, one which includes the poles, it calculates the winds on a grid offset by half a grid point in longitude and latitude, test 3 Starting with analytically generated winds on an evenly spaced grid, it calculates the winds on a gaussian grid. For each test an rms error is calculated and written to the screen. In addition, the output to the screen is placed in the file screen.asc. The numerical results are written to a series of Netcdf files. The contents of these files are described in the screen output. After completing the tests, documentation is written to the file spheremodule.doc. DOCUMENTATION Documentation written to the file spheremodule.doc can be obtained without running the tests after importing the spheretest module by typing spheretest.document() A brief view of the documentation consisting of the overview can be written to the file spheremodule.doc after importing the sphere module by typing spheretest.document(brief = 'yes') """ import sys, string, sphere import spherepack, numpy, math #spherepack.set_pyfort_option(spherepack.MIRROR) debug = 0 # set to 1 for debug prints radius = 6.37122e06 writeTestcase = 'yes' try: import cdms2 except ImportError: print 'Can not write test case results to netCDF files without module cdms2' writeTestcase = 'no' def document(brief = 'no'): #---------------------------------------------------------------------------------------- # # purpose: 'document' writes the doc strings contained in the sphere module # to a file as documentation for the user # # usage: import sphere # sphere.document() # # passed : nothing # # returned: nothing # #---------------------------------------------------------------------------------------- std = sys.stdout sys.stdout = open('spheremodule.doc', 'w') print '**********************************************************************************************\n' print '****************** Overview of the CDAT interface to the NCAR SPHEREPACK 3.0 *****************\n' print '**********************************************************************************************\n' print sphere.__doc__ print print if brief != 'no': return None print ' ************************************************************************************\n' print ' ************************** Sphere class documentation ******************************\n' print ' ************************************************************************************\n' sphere.help('Sphere') # how to make an instance of a class print print ' **************************** Sphere class functions ********************************\n' SphereList = ['div', 'idiv', 'vrt', 'ivrt', 'idvt', 'vts', 'grad', 'igrad', 'slap', 'islap', 'vlap', 'ivlap', 'sfvp', 'isfvp', 'truncation', 'sha', 'shs', 'vha', 'vhs'] for name in SphereList: command = "sphere.help(" + "name" + ")" exec command print print ' ************************************************************************************\n' print ' **************************** Regrid class documentation ****************************\n' print ' ************************************************************************************\n' sphere.help('Regrid') # how to make an instance of a class print print ' ****************************** Regrid class functions ******************************\n' sphere.help('regridScalar') print sphere.help('regridVector') print print ' ************************************************************************************\n' print ' *************************** Shiftgrid class documentation **************************\n' print ' ************************************************************************************\n' sphere.help('Shiftgrid') # how to make an instance of a class print print ' ***************************** Shiftgrid class functions ****************************\n' sphere.help('shiftScalar') print sphere.help('shiftVector') print print ' ************************************************************************************\n' print ' **************************** Utility documentation *********************************\n' print ' ************************************************************************************\n' sphere.help('gridGenerator') # how to make an instance of a class print sphere.help('truncate') print sys.stdout = std return None # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # ++++++++++++++++++++++++++++++++ Autotest Functions +++++++++++++++++++++++++++++++++ # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def sfvp(): #----------------------------------------------------------------------------------- # # purpose: starting with analytically generated winds, calculate the stream # function and velocity potential using Spherepack and compare with # analytically generated stream function and velocity potential. # # usage: sfvp() # # passed : nothing # # returned: nothing # #----------------------------------------------------------------------------------- sendmsg('**************** calculate the stream function and velocity potential on gaussian grid *****************') sendmsg(' ') testError = 0 comp = 'computed' nlon = 128 nlat = 64 lonvals, latvals, timevals, u, v, sfexact, vpexact = sphere_test(nlon, nlat, 'v', 'gaussian') nt = len(timevals) x = sphere.Sphere(lonvals, latvals, numberLevels = 0, numberTimes = nt, computed_stored = comp) sfcal, vpcal = x.sfvp(u, v) scale = radius # scale exact functions to radius for the earth sfexact = sphere.geoscale(scale, sfexact) vpexact = sphere.geoscale(scale, vpexact) sfexact = remove_offset(sfcal, sfexact) # subtract the offset vpexact = remove_offset(vpcal, vpexact) sendmsg('******* compare results') rms = rmserror(sfcal, sfexact) # stream function rms error sendmsg('expected normalized rms error in stream function computation is less than 1.e-05') sendmsg('calculated normalized rms error in stream function computation =', rms) sendmsg(' ') if rms > 1.e-05: testError = testError + 1 rms = rmserror(vpcal, vpexact) # velocity potential rms error sendmsg( 'expected normalized rms error in velocity potential computation is less than 1.e-05') sendmsg( 'calculated normalized rms error in velocity potential computation =', rms) sendmsg(' ') if rms > 1.e-05: testError = testError + 1 if writeTestcase == 'yes': sendmsg('******* write data') sendmsg( 'calculated stream function written to sfcal.nc') # write netcdf file writeField(lonvals, latvals, timevals, 'sfcal', sfcal) sendmsg( 'calculated velocity potential written to vpcal.nc') writeField(lonvals, latvals, timevals, 'vpcal', vpcal) sendmsg(' ') sendmsg( 'exact stream function written to sfexact.nc') # write netcdf file writeField(lonvals, latvals, timevals, 'sfexact', sfexact) sendmsg( 'exact velocity potential written to vpexact.nc') writeField(lonvals, latvals, timevals, 'vpexact', vpexact) sendmsg(' ') sendmsg(' ') return testError def shift(): #----------------------------------------------------------------------------------- # # purpose: starting with analytically generate winds on a regular grid, one # which includes the poles, calculate the winds on a grid offset by # half a grid point in longitude and latitude, # # usage: shift() # # passed : nothing # # returned: nothing # #----------------------------------------------------------------------------------- sendmsg('************ vector shift from a regular evenly spaced grid and a offset evenly spaced grid ************') sendmsg(' ') testError = 0 comp = 'computed' nlon = 144 nlat = 72 lonvals, latvals, timevals, u, v = vshift_testfunction(nlon, nlat) offlonvals, offlatvals, timevals, uexact, vexact = vshift_testfunction(nlon, nlat, grid_type = 'offset') nt = len(timevals) x = sphere.Shiftgrid(lonvals, latvals, numberTimes = nt) ucal, vcal = x.shiftVector(u,v) sendmsg('******* compare results') rms = rmserror(ucal, uexact) # zonal wind rms error sendmsg( 'expected normalized error in zonal wind computation is less than 1.e-06') sendmsg( 'calculated normalized rms error in zonal wind computation =', rms) sendmsg( ' ') if rms > 1.e-06: testError = testError + 1 rms = rmserror(vcal, vexact) sendmsg( 'expected normalized rms error in meridional wind computation is less than 1.e-06') sendmsg( 'calculated normalized rms error in meridional wind computation =', rms) sendmsg(' ') if rms > 1.e-06: testError = testError + 1 if writeTestcase == 'yes': sendmsg( '******* write data') sendmsg( 'calculated zonal wind on the offset grid written to uoffcal.nc') # write netcdf file writeField(offlonvals, offlatvals, timevals, 'uoffcal', ucal) sendmsg( 'calculated meridional wind on the offset grid wriiten to voffcal.nc') writeField(offlonvals, offlatvals, timevals, 'voffcal', vcal) sendmsg(' ') sendmsg( 'exact zonal wind on the offset grid written to uoffexact.nc') # write netcdf file writeField(offlonvals, offlatvals, timevals, 'uoffexact', uexact) sendmsg( 'exact meritional wind on the offset grid written to voffexact.nc') writeField(offlonvals, offlatvals, timevals, 'voffexact', vexact) sendmsg(' ') sendmsg(' ') return testError def regrid(): #----------------------------------------------------------------------------------- # # purpose: starting with analytically generate winds on an evenly spaced grid, # calculate the winds on a gaussian grid. # # usage: regrid() # # passed : nothing # # returned: nothing # #----------------------------------------------------------------------------------- sendmsg('******************************* vector regridding on a sphere ***********************************') sendmsg(' ') testError = 0 nlone = 36 nlate = 19 lonArrayIn, latArrayIn, timevals, uIn, vIn = regrid_testfunction(nlone, nlate, 'v', 'even') nlong = 128 nlatg = 64 lonArrayOut, latArrayOut, timevals, uOut, vOut = regrid_testfunction(nlong, nlatg, 'v', 'gaussian') nt = len(timevals) x = sphere.Regrid(lonArrayOut, latArrayOut, lonArrayIn, latArrayIn, numberTimes = nt) ucal, vcal = x.regridVector(uIn, vIn) sendmsg('******* compare results') rms = rmserror(uOut, ucal) sendmsg( 'expected normalized rms error in regridded zonal wind is less than 1.e-06') sendmsg( 'calculated normalized rms error in regridded zonal wind =', rms) sendmsg(' ') if rms > 1.e-06: testError = testError + 1 rms = rmserror(vOut, vcal) sendmsg( 'expected normalized rms error in regridded meridonal wind is less than 1.e-06') sendmsg( 'calculted normalized rms error in regridded meridonal wind =', rms) sendmsg(' ') if rms > 1.e-06: testError = testError + 1 if writeTestcase == 'yes': sendmsg('******* write data') sendmsg( 'calculated zonal wind on the new grid written to uregridcal.nc') # write netcdf file writeField(lonArrayOut, latArrayOut, timevals, 'uregridcal', ucal) sendmsg( 'calculated meridional wind on the new grid written to vregridcal.nc') writeField(lonArrayOut, latArrayOut, timevals, 'vregridcal', vcal) sendmsg(' ') sendmsg( 'exact zonal wind on the new grid written to uregridexact.nc') # write netcdf file writeField(lonArrayOut, latArrayOut, timevals, 'uregridexact', uOut) sendmsg( 'exact meridional wind on the new grid written to vregridexact.nc') writeField(lonArrayOut, latArrayOut, timevals, 'vregridexact', vOut) sendmsg(' ') sendmsg( 'original zonal wind on the grid written to uorig.nc') writeField(lonArrayIn, latArrayIn, timevals, 'uorig', uIn) sendmsg( 'original meridional wind on the grid written to vorig.nc') writeField(lonArrayIn, latArrayIn, timevals, 'vorig', vIn) return testError def writeField(lons, lats, tmes, varname, dataField): #----------------------------------------------------------------------------------- # # purpose: write an output field # # usage: # # passed : lons, lats, tmes, filename, varname, dataField # # returned: None # #----------------------------------------------------------------------------------- fileObj = cdms2.createDataset(varname + '.nc') lon_axis = fileObj.createAxis('longitude', lons) lon_axis.units = "degrees_east" lat_axis = fileObj.createAxis('latitude', lats) lat_axis.units = "degrees_north" tme_axis = fileObj.createAxis('time', tmes) tme_axis.units = "months" var = fileObj.createVariable(varname, numpy.float32, (tme_axis, lat_axis, lon_axis)) # variable without data var[:] = dataField # copy in the data fileObj.close() return None def writeasc(name, a): #----------------------------------------------------------------------------------- # # purpose: write the data into an ascii file # # passed : name - filename prefix # a - data # # returned: return # #----------------------------------------------------------------------------------- r = numpy.ravel(a) output = open(name + '.asc', 'w') Format = '%12.4E' count = 0 for item in r: output.write(Format % (item,)) count = count + 1 if count == 8: output.write('\n') count = 0 output.close() return None def sendmsg(msg, value = None, screen = 'no'): #------------------------------------------------------------------------------ # # purpose: send the same message to the screen and to a file # # passed : msg - the string # # returned: return # #------------------------------------------------------------------------------ if value is None: if screen != 'no': print msg output.write(msg + '\n') else: if screen != 'no': print msg, `value` output.write(msg + ' %15.11e\n' % (value,)) return None def remove_offset(cal, test): #----------------------------------------------------------------------------------- # # purpose: shift the array so that the test array has an approximate zero mean to agree # with the calculated one. # # passed : the two data sets # # returned: test # #----------------------------------------------------------------------------------- dif = numpy.ravel(cal) - numpy.ravel(test) # remove offset between the test and calculated arrays offset = numpy.sum(dif)/len(dif) test = test + offset test = test.astype(numpy.float32) return test def rmserror(data1, data2): #----------------------------------------------------------------------------------- # # purpose: compute the rms error for two data sets having the same shape # # passed : the two data sets # # returned: rms error # #----------------------------------------------------------------------------------- if data1.shape != data2.shape: print 'Error in shape in rmserror' raise ValueError d1 = numpy.ravel(data1) d2 = numpy.ravel(data2) sq = d1*d1 # find average magnitude avg = numpy.sqrt( numpy.sum(sq)/len(d1)) sq = (d1 - d2)*(d1 - d2) error = numpy.sum(sq)/len(d1) rmserror = (numpy.sqrt(error))/avg return rmserror def sphere_test(nlon, nlat, sorv = 's', grid_choice = 'even'): #----------------------------------------------------------------------------------- # # purpose: compute the scalar and vector test functions used in testrssph.f # and testtrvsph.f in geophysical coordinates # #----------------------------------------------------------------------------------- if sorv != 's' and sorv != 'v': print 'Must pass s or v to select scalar or vector' raise ValueError return # generate the correct geophysical grid points lonlistrad = [] # longitudes in radians delta = 2.*math.pi/nlon for i in range(nlon): value = i*delta lonlistrad.append(value) lons = numpy.array(lonlistrad, numpy.float64) # in radians for computation lonlist = map( (lambda x: (180./math.pi)*x), lonlistrad) lonvals = numpy.array(lonlist, numpy.float64) # in degrees for return if grid_choice == 'even': latlistrad = [] # latitudes in radians delta = math.pi/(nlat - 1) for i in range(nlat): value = math.pi/2. - i*delta latlistrad.append(value) lats = numpy.array(latlistrad, numpy.float64) # in radians for computation latlist = map( (lambda x: (180./math.pi)*x), latlistrad) latvals = numpy.array(latlist, numpy.float64) # array of latitudes in degrees elif grid_choice == 'gaussian': ldwork = nlat*(nlat + 2) work = numpy.zeros((ldwork,),'d') points, wts, ierror = spherepack.gaqd(nlat, work) # get colatitudes from gaqd if ierror != 0: print 'In return from call to gaqd ierror = ', ierror raise ValueError # convert points to latitudes colatlistrad = list(points) latlistrad = map( (lambda x: math.pi/2. - x), colatlistrad) # convert to latitudes lats = numpy.array(latlistrad, numpy.float64) # in radians for computation latlist = map( (lambda x: (180./math.pi)*x), latlistrad) latvals = numpy.array(latlist, numpy.float64) else: print 'grid_choice must be even or gaussian' raise ValueError return # generate the data on geophysical grid points timevals = numpy.array( [0.0], numpy.float) if sorv == 's': sf = numpy.zeros((1,nlat,nlon), numpy.float) # malloc in c order for i in range(len(lons)): # calculate scalar test function p = lons[i] cosp = math.cos(p) sinp = math.sin(p) for j in range(len(lats)): t = lats[j] cost = math.cos(t) sint = math.sin(t) x = cost*cosp y = cost*sinp z = sint sf[0,j,i] = math.exp(x*y*z) # c order sf = sf.astype(numpy.float32) return lonvals, latvals, timevals, sf else: u = numpy.zeros((1,nlat,nlon), numpy.float) # malloc in c order v = numpy.zeros((1,nlat,nlon), numpy.float) sf = numpy.zeros((1,nlat,nlon), numpy.float) vp = numpy.zeros((1,nlat,nlon), numpy.float) for i in range(len(lons)): # calculate vector test functions p = lons[i] cosp = math.cos(p) sinp = math.sin(p) for j in range(len(lats)): t = lats[j] cost = math.cos(t) sint = math.sin(t) x = cost*cosp y = cost*sinp z = sint ex = math.exp(x) ey = math.exp(y) ez = math.exp(z) emz = math.exp(-z) u[0,j,i] = -ex*sinp + ey*sint*sinp + emz*cost # c order v[0,j,i] = -( ex*sint*cosp - ey*cosp - ez*cost ) sf[0,j,i] = ey + emz vp[0,j,i] = ex + ez u = u.astype(numpy.float32) v = v.astype(numpy.float32) sf = sf.astype(numpy.float32) vp = vp.astype(numpy.float32) return lonvals, latvals, timevals, u, v, sf, vp def regrid_testfunction(nlon, nlat, sorv = 'h', grid_choice = 'even'): #----------------------------------------------------------------------------------- # # purpose: compute the scalar and vector test functions used in testrssph.f # and testtrvsph.f in geophysical coordinates # # usage: # # #----------------------------------------------------------------------------------- if sorv != 's' and sorv != 'v': print 'Must pass s or v to select scalar or vector' raise ValueError return # generate the correct geophysical grid points lonlistrad = [] # longitudes in radians delta = 2.*math.pi/nlon for i in range(nlon): value = i*delta lonlistrad.append(value) lons = numpy.array(lonlistrad, numpy.float64) # in radians for computation lonlist = map( (lambda x: (180./math.pi)*x), lonlistrad) lonvals = numpy.array(lonlist, numpy.float64) # in degrees for return if grid_choice == 'even': latlistrad = [] # latitudes in radians delta = math.pi/(nlat - 1) for i in range(nlat): value = math.pi/2. - i*delta latlistrad.append(value) lats = numpy.array(latlistrad, numpy.float64) # in radians for computation latlist = map( (lambda x: (180./math.pi)*x), latlistrad) latvals = numpy.array(latlist, numpy.float64) # array of latitudes in degrees elif grid_choice == 'gaussian': ldwork = nlat*(nlat + 2) work = numpy.zeros((ldwork,),'d') points, wts, ierror = spherepack.gaqd(nlat, work) # get colatitudes from gaqd if ierror != 0: print 'In return from call to gaqd ierror = ', ierror raise ValueError # convert points to latitudes colatlistrad = list(points) latlistrad = map( (lambda x: math.pi/2. - x), colatlistrad) # convert to latitudes lats = numpy.array(latlistrad, numpy.float64) # in radians for computation latlist = map( (lambda x: (180./math.pi)*x), latlistrad) latvals = numpy.array(latlist, numpy.float64) else: print 'grid_choice must be even or gaussian' raise ValueError return # generate the data on geophysical grid points timevals = numpy.array( [0.0], numpy.float) if sorv == 's': sf = numpy.zeros((1,nlat,nlon), numpy.float) # malloc in c order for i in range(len(lons)): # calculate scalar test function p = lons[i] cosp = math.cos(p) sinp = math.sin(p) for j in range(len(lats)): t = lats[j] cost = math.cos(t) sint = math.sin(t) x = cost*cosp y = cost*sinp z = sint sf[0,j,i] = math.exp(x*y*z) # c order sf = sf.astype(numpy.float32) return lonvals, latvals, timevals, sf else: u = numpy.zeros((1,nlat,nlon), numpy.float) # malloc in c order v = numpy.zeros((1,nlat,nlon), numpy.float) for i in range(len(lons)): # calculate vector test functions p = lons[i] cosp = math.cos(p) sinp = math.sin(p) for j in range(len(lats)): t = lats[j] cost = math.cos(t) sint = math.sin(t) x = cost*cosp y = cost*sinp z = sint ex = math.exp(x) ey = math.exp(y) ez = math.exp(z) emz = math.exp(-z) u[0,j,i] = -ex*sinp + ey*sint*sinp + emz*cost # c order v[0,j,i] = -( ex*sint*cosp - ey*cosp - ez*cost ) u = u.astype(numpy.float32) v = v.astype(numpy.float32) return lonvals, latvals, timevals, u, v def sshift_testfunction(nlon, nlat, grid_type = 'regular'): #----------------------------------------------------------------------------------- # # purpose: compute the scalar test function used in testsshifte.f # in geophysical coordinates # # usage: # # #----------------------------------------------------------------------------------- if grid_type != 'regular' and grid_type != 'offset': print 'Must pass regular or offset to select grid from sshift_testfunction' raise ValueError return # ------- generate the longitude grid points reg_lonlistrad = [] # regular longitudes in radians off_lonlistrad = [] # offset longitudes in radians delta = 2.*math.pi/nlon shiftlon = delta/2. for i in range(nlon): value = i*delta svalue = shiftlon + value reg_lonlistrad.append(value) off_lonlistrad.append(svalue) reglons = numpy.array(reg_lonlistrad, numpy.float64) # for use in computation offlons = numpy.array(off_lonlistrad, numpy.float64) # for use in computation lonlist = map( (lambda x: (180./math.pi)*x), reg_lonlistrad) # degrees for return reglonvals = numpy.array(lonlist, numpy.float64) lonlist = map( (lambda x: (180./math.pi)*x), off_lonlistrad) offlonvals = numpy.array(lonlist, numpy.float64) # generate the data on geophysical grid points timevals = numpy.array( [0.0], numpy.float) if grid_type == 'regular': # -------- generate the regular latitude grid points latlistrad = [] # latitudes in radians delta = math.pi/nlat for j in range(nlat + 1): value = math.pi/2. - j*delta latlistrad.append(value) lats = numpy.array(latlistrad, numpy.float64) # latitudes for computation latlist = map( (lambda x: (180./math.pi)*x), latlistrad) latvals = numpy.array(latlist, numpy.float64) # array of latitudes in degrees sf = numpy.zeros((1,nlat + 1,nlon), numpy.float32) # malloc in c order for i in range(len(reglons)): # calculate scalar test function p = reglons[i] cosp = math.cos(p) sinp = math.sin(p) for j in range(len(lats)): t = lats[j] cost = math.cos(t) sint = math.sin(t) x = cost*cosp y = cost*sinp z = sint sf[0,j,i] = math.exp(x + y + z) # c order sf = sf.astype(numpy.float32) return reglonvals, latvals, timevals, sf else: # -------- generate the offset latitude grid points latlistrad = [] # latitudes in radians delta = math.pi/nlat shiftlat = delta/2. for j in range(nlat): value = math.pi/2. - shiftlat - j*delta latlistrad.append(value) lats = numpy.array(latlistrad, numpy.float64) # latitudes for computation latlist = map( (lambda x: (180./math.pi)*x), latlistrad) latvals = numpy.array(latlist, numpy.float64) # array of latitudes in degrees sf = numpy.zeros((1,nlat,nlon), numpy.float32) # malloc in c order for i in range(len(offlons)): # calculate scalar test function p = offlons[i] cosp = math.cos(p) sinp = math.sin(p) for j in range(len(lats)): t = lats[j] cost = math.cos(t) sint = math.sin(t) x = cost*cosp y = cost*sinp z = sint sf[0,j,i] = math.exp(x + y + z) # c order sf = sf.astype(numpy.float32) return offlonvals, latvals, timevals, sf def vshift_testfunction(nlon, nlat, grid_type = 'regular'): #----------------------------------------------------------------------------------- # # purpose: compute the vector test function used in testvshifte.f # in geophysical coordinates # # usage: # # #----------------------------------------------------------------------------------- if grid_type != 'regular' and grid_type != 'offset': print 'Must pass regular or offset to select grid from vshift_testfunction' raise ValueError return # ------- generate the longitude grid points reg_lonlistrad = [] # regular longitudes in radians off_lonlistrad = [] # offset longitudes in radians delta = 2.*math.pi/nlon shiftlon = delta/2. for i in range(nlon): value = i*delta svalue = shiftlon + value reg_lonlistrad.append(value) off_lonlistrad.append(svalue) reglons = numpy.array(reg_lonlistrad, numpy.float64) # for use in computation offlons = numpy.array(off_lonlistrad, numpy.float64) # for use in computation lonlist = map( (lambda x: (180./math.pi)*x), reg_lonlistrad) # degrees for return reglonvals = numpy.array(lonlist, numpy.float64) lonlist = map( (lambda x: (180./math.pi)*x), off_lonlistrad) offlonvals = numpy.array(lonlist, numpy.float64) # generate the data on geophysical grid points timevals = numpy.array( [0.0], numpy.float) if grid_type == 'regular': # -------- generate the regular latitude grid points latlistrad = [] # latitudes in radians delta = math.pi/nlat for j in range(nlat + 1): value = math.pi/2. - j*delta latlistrad.append(value) lats = numpy.array(latlistrad, numpy.float64) # latitudes for computation latlist = map( (lambda x: (180./math.pi)*x), latlistrad) latvals = numpy.array(latlist, numpy.float64) # array of latitudes in degrees u = numpy.zeros((1,nlat + 1,nlon), numpy.float32) # malloc in c order v = numpy.zeros((1,nlat + 1,nlon), numpy.float32) for i in range(len(reglons)): # calculate scalar test function p = reglons[i] cosp = math.cos(p) sinp = math.sin(p) for j in range(len(lats)): t = lats[j] cost = math.cos(t) sint = math.sin(t) x = cost*cosp y = cost*sinp z = sint ex = math.exp(x) ey = math.exp(y) ez = math.exp(z) emz = math.exp(-z) u[0,j,i] = -ex*sinp + ey*sint*sinp + emz*cost # c order v[0,j,i] = -ex*sint*cosp + ey*cosp + ez*cost u = u.astype(numpy.float32) v = v.astype(numpy.float32) return reglonvals, latvals, timevals, u, v else: # -------- generate the offset latitude grid points latlistrad = [] # latitudes in radians delta = math.pi/nlat shiftlat = delta/2. for j in range(nlat): value = math.pi/2. - shiftlat - j*delta latlistrad.append(value) lats = numpy.array(latlistrad, numpy.float64) # latitudes for computation latlist = map( (lambda x: (180./math.pi)*x), latlistrad) latvals = numpy.array(latlist, numpy.float64) # array of latitudes in degrees u = numpy.zeros((1,nlat,nlon), numpy.float32) # malloc in c order v = numpy.zeros((1,nlat,nlon), numpy.float32) for i in range(len(offlons)): # calculate scalar test function p = offlons[i] cosp = math.cos(p) sinp = math.sin(p) for j in range(len(lats)): t = lats[j] cost = math.cos(t) sint = math.sin(t) x = cost*cosp y = cost*sinp z = sint ex = math.exp(x) ey = math.exp(y) ez = math.exp(z) emz = math.exp(-z) u[0,j,i] = -ex*sinp + ey*sint*sinp + emz*cost # c order v[0,j,i] = -ex*sint*cosp + ey*cosp + ez*cost u = u.astype(numpy.float32) v = v.astype(numpy.float32) return offlonvals, latvals, timevals, u, v if __name__ == "__main__": output = open('test.asc', 'w') # global file name print 'Running the test computations' te1 = sfvp() te2 = shift() te3 = regrid() testError = te1 + te2 + te3 write = document() sendmsg(' ') sendmsg('*********') sendmsg('General information on the use of SPHERPACK has been written to the file spheremodule.doc.') sendmsg('*********') sendmsg(' ') if testError == 0: print 'Testing Completed Successfully' else: print 'Testing Completed But It May Have Problems' print 'Some details on the testing have been written to the file test.asc.' print 'General information on the use of SPHEREPACK has been written to the file spheremodule.doc.' output.close()