[cig-commits] r7650 - in mc/1D/hc/trunk: . hcplates
becker at geodynamics.org
becker at geodynamics.org
Thu Jul 12 12:08:53 PDT 2007
Author: becker
Date: 2007-07-12 12:08:51 -0700 (Thu, 12 Jul 2007)
New Revision: 7650
Modified:
mc/1D/hc/trunk/.gmtcommands
mc/1D/hc/trunk/Makefile
mc/1D/hc/trunk/README.TXT
mc/1D/hc/trunk/expansion.par
mc/1D/hc/trunk/ggrd.h
mc/1D/hc/trunk/ggrd_grdtrack_util.c
mc/1D/hc/trunk/ggrd_grdtrack_util.h
mc/1D/hc/trunk/ggrd_readgrds.c
mc/1D/hc/trunk/ggrd_struc.h
mc/1D/hc/trunk/ggrd_velinterpol.c
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc2gmt
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/hc_input.c
mc/1D/hc/trunk/hc_misc.c
mc/1D/hc/trunk/hc_polsol.c
mc/1D/hc/trunk/hc_solve.c
mc/1D/hc/trunk/hcplates/Makefile
mc/1D/hc/trunk/hcplates/findplate.c
mc/1D/hc/trunk/hcplates/hc_ptrot_init.c
mc/1D/hc/trunk/hcplates/hcplates_init.c
mc/1D/hc/trunk/hcplates/ptrot.c
mc/1D/hc/trunk/hcplates/run_example.mk
mc/1D/hc/trunk/main.c
mc/1D/hc/trunk/pow.sh_ana
mc/1D/hc/trunk/pow.shana
mc/1D/hc/trunk/rick_sh_c.c
mc/1D/hc/trunk/sh.h
mc/1D/hc/trunk/sh_ana.c
mc/1D/hc/trunk/sh_exp.c
mc/1D/hc/trunk/sh_model.c
mc/1D/hc/trunk/sh_power.c
mc/1D/hc/trunk/sh_rick.h
mc/1D/hc/trunk/sh_shana.h
mc/1D/hc/trunk/sh_syn.c
mc/1D/hc/trunk/shana_sh.c
mc/1D/hc/trunk/test_sh_routines
Log:
Minor updates in the hc part, some modifications to the hcplates compilation scheme.
Modified: mc/1D/hc/trunk/.gmtcommands
===================================================================
--- mc/1D/hc/trunk/.gmtcommands 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/.gmtcommands 2007-07-12 19:08:51 UTC (rev 7650)
@@ -1,3 +1,3 @@
# GMT common arguments shelf
--R000.000000/359.000000/-89.500000/089.500000
+-R0/360/-90/90
EOF
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/Makefile 2007-07-12 19:08:51 UTC (rev 7650)
@@ -1,11 +1,15 @@
#
#
# makefile for experimental Hager & O'Connell routines
+# and hcplates Ricard/Vigny type plate velocity inversions
#
+# see source files for comments and reference to original authors
+#
+#
#
#
LD = $(CC)
-CFLAGS = $(CFLAGS_DEBUG) -Wall
+CFLAGS = $(CFLAGS_DEBUG)
#
# EDIT HERE FOR GMT VERSION
@@ -135,17 +139,21 @@
all: dirs libs hc hc.dbg hc_extract_sh_layer \
sh_syn sh_ana sh_power \
- ggrd_test
+ ggrd_test hcplates
libs: dirs hc_lib $(HEAL_LIBS) $(RICK_LIB)
hc_lib: $(HC_LIBS) $(GGRD_LIBS)
-really_all: dirs proto all
+really_all: proto all
proto: hc_auto_proto.h
+hcplates:
+ cd hcplates; \
+ make ;\
+ cd ..
sh_test: $(LIBS) $(INCS) $(ODIR)/sh_test.o
Modified: mc/1D/hc/trunk/README.TXT
===================================================================
--- mc/1D/hc/trunk/README.TXT 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/README.TXT 2007-07-12 19:08:51 UTC (rev 7650)
@@ -4,6 +4,13 @@
tools and libraries installed. See the Makefile for comments and what
to modify.
+
+What is described below is the hc Hager & O'Connell (1981) forward
+flow and geoid computation. For plate velocity inversions following
+the Ricard & Vigny (1989) method, see the hcplates subdirectory and
+the README there.
+
+
The Makefile assumes that the following environment variables are
predefined:
@@ -88,6 +95,10 @@
SPHERICAL HARMONICS FORMAT
+
+(A) Regular (long) format, which allows for both scalar and vector
+harmonics
+
All single layer spherical harmonics are in the fully normalized,
physical convention, e.g. Dahlen & Tromp (1998).
@@ -132,7 +143,14 @@
...
+(B) Short format
+As above, but will only expect a single integer, lmax, as the header
+for a spherical harmonic expansion.
+
+
+
+
OTHER INPUT FILE FORMATS
1) Viscosity structure (hc_assign_viscosity)
@@ -162,35 +180,46 @@
a) extract SH expansion of radial velocity at layer 2
- hc_extract_sh_layer sol.bin 2 1
+ hc_extract_sh_layer vel.sol.bin 2 1
b) extract SH expansion of radial velocity at layer 2 and convert
to spatial expansion
- hc_extract_sh_layer sol.bin 2 1 0 | sh_syn
+ hc_extract_sh_layer vel.sol.bin 2 1 0 | sh_syn
c) extract SH expansion of poloidal and toroidal velocity at layer 5 and convert
to spatial expansion of v_theta v_phi
- hc_extract_sh_layer sol.bin 5 2 0 | sh_syn
+ hc_extract_sh_layer vel.sol.bin 5 2 0 | sh_syn
+USING THE SPHERICAL HARMONICS TOOLS AS STANDALONE
-DIRECTORY LISTING:
+1) Convert scalar values from a GMT/Netcdf grd file into a spherical
+ harmonics expansion of degree 31
-becker at jackie:~/progs/src/hc-svn > ls *.c *.h
-ggrd_grdtrack_util.c hc.h prem.h sh_rick.h
-ggrd_grdtrack_util.h hc_init.c prem_util.c sh_shana.h
-ggrd.h hc_input.c rick_fft_c.c sh_spherepack.h
-ggrd_readgrds.c hc_matrix.c rick_sh_c.c sh_syn.c
-ggrd_struc.h hc_misc.c sh_ana.c sh_test.c
-ggrd_test.c hc_output.c shana_sh.c simple_test.c
-ggrd_velinterpol.c hc_polsol.c sh_exp.c spherepack_sh.c
-hc_auto_proto.h hc_propagator.c sh.h test_fft.c
-hc_auto_proto.temp.h hc_solve.c sh_model.c
-hc_extract_sh_layer.c hc_torsol.c sh_power.c
-hc_filenames.h main.c sh_rick_ftrn.h
+ a) obtain scalara data at the Gaussian intergration points
+
+ sh_ana -127 | grdtrack -Lx -G$datadir/etopo5/etopo5.1.grd > etopo5.dat
+
+
+ b) expand
+
+ cat etopo5.dat | sh_ana 127 > etopo5.ab
+
+2) convert spherical harmonics to spatial expansion
+
+ cat etopo5.ab | sh_syn > etopo5.127.dat
+
+
+
+
+
+
+
+
+
Modified: mc/1D/hc/trunk/expansion.par
===================================================================
--- mc/1D/hc/trunk/expansion.par 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/expansion.par 2007-07-12 19:08:51 UTC (rev 7650)
@@ -1 +1 @@
- 0.0000000000 1.0000000000 360.0000000000 -90.0000000000 1.0000000000 90.0000000000 361 181
+ 0.0000000000 2.0000000000 360.0000000000 -90.0000000000 2.0000000000 90.0000000000 181 91
Modified: mc/1D/hc/trunk/ggrd.h
===================================================================
--- mc/1D/hc/trunk/ggrd.h 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/ggrd.h 2007-07-12 19:08:51 UTC (rev 7650)
@@ -13,6 +13,7 @@
*/
#define GGRD_CPREC double
#define GGRD_EPS 5e-15
+#define ggrd_boolean unsigned char
#define HC_FLT_FORMAT "%lf"
#endif
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -35,39 +35,39 @@
grdfile: filename for 2-D and prefix for 3-D
depth_file: file with depth layers for 3-D
edge_info_string: as in GMT
- g: ggrd_gt structure
+ g: ggrd_gt structure, pass allocated
returns error code
*/
-int ggrd_grdtrack_init_general(unsigned char is_three,
+int ggrd_grdtrack_init_general(ggrd_boolean is_three,
char *grdfile, char *depth_file,
char *gmt_edgeinfo_string,
struct ggrd_gt *g,
- unsigned char verbose,
- unsigned char change_z_sign)
+ ggrd_boolean verbose,
+ ggrd_boolean change_z_sign)
{
- static unsigned char bilinear = FALSE; /* cubic by default */
- static double west=0.,east=0.,south=0.,north=0.; /* whole region by default */
+ /* this is a switch and can be left in */
+ static ggrd_boolean bilinear = TRUE; /* bilinear by default */
+
int pad[4];
- unsigned char geographic_in; /* this is set by grdtrack_init */
+ ggrd_boolean geographic_in; /* this is set by grdtrack_init */
int i,j;
float zavg;
- if(g->init){
- fprintf(stderr,"ggrd_grdtrack_init_general: error: call only once\n");
- return 1;
- }
+ /* clear all entries */
+ g->east = g->west = g->south = g->north = 0.0;
g->is_three = is_three;
+ g->init = FALSE;
#ifdef USE_GMT4
- if(ggrd_grdtrack_init(&west,&east,&south,&north,&g->f,&g->mm,
+ if(ggrd_grdtrack_init(&g->west,&g->east,&g->south,&g->north,&g->f,&g->mm,
grdfile,&g->grd,&g->edgeinfo,
gmt_edgeinfo_string,&geographic_in,
pad,is_three,depth_file,&g->z,&g->nz,
bilinear,verbose,change_z_sign,&g->bcr))
return 2;
#else
- if(ggrd_grdtrack_init(&west,&east,&south,&north,&g->f,&g->mm,
+ if(ggrd_grdtrack_init(&g->west,&g->east,&g->south,&g->north,&g->f,&g->mm,
grdfile,&g->grd,&g->edgeinfo,
gmt_edgeinfo_string,&geographic_in,
pad,is_three,depth_file,&g->z,&g->nz,
@@ -128,9 +128,9 @@
*/
int ggrd_grdtrack_rescale(struct ggrd_gt *g,
- unsigned char take_log10, /* take log10() */
- unsigned char take_power10, /* take 10^() */
- unsigned char rescale, /* rescale? */
+ ggrd_boolean take_log10, /* take log10() */
+ ggrd_boolean take_power10, /* take 10^() */
+ ggrd_boolean rescale, /* rescale? */
double scale /* factor for rescaling */)
{
int i,j,k;
@@ -162,12 +162,12 @@
interpolation wrapper, uses r, theta, phi input. return value and TRUE if success,
undefined and FALSE else
*/
-unsigned char ggrd_grdtrack_interpolate_rtp(double r,double t,double p,
+ggrd_boolean ggrd_grdtrack_interpolate_rtp(double r,double t,double p,
struct ggrd_gt *g,double *value,
- unsigned char verbose)
+ ggrd_boolean verbose)
{
double x[3];
- unsigned char result;
+ ggrd_boolean result;
if(!g->init){ /* this will not necessarily work */
fprintf(stderr,"ggrd_grdtrack_interpolate_rtp: error, g structure not initialized\n");
return FALSE;
@@ -205,12 +205,12 @@
undefined and FALSE else
this mean lon lat z
*/
-unsigned char ggrd_grdtrack_interpolate_xyz(double x,double y,double z,
+ggrd_boolean ggrd_grdtrack_interpolate_xyz(double x,double y,double z,
struct ggrd_gt *g,double *value,
- unsigned char verbose)
+ ggrd_boolean verbose)
{
double xloc[3];
- unsigned char result;
+ ggrd_boolean result;
if(!g->init){ /* this will not necessarily work */
fprintf(stderr,"ggrd_grdtrack_interpolate_xyz: error, g structure not initialized\n");
return FALSE;
@@ -247,13 +247,13 @@
undefined and FALSE else
*/
-unsigned char ggrd_grdtrack_interpolate_tp(double t,double p,
+ggrd_boolean ggrd_grdtrack_interpolate_tp(double t,double p,
struct ggrd_gt *g,
double *value,
- unsigned char verbose)
+ ggrd_boolean verbose)
{
double x[3];
- unsigned char result;
+ ggrd_boolean result;
if(!g->init){ /* this will not necessarily work */
fprintf(stderr,"ggrd_grdtrack_interpolate_tp: error, g structure not initialized\n");
return FALSE;
@@ -289,13 +289,13 @@
return value and TRUE if success,
undefined and FALSE else
*/
-unsigned char ggrd_grdtrack_interpolate_xy(double xin,double yin,
+ggrd_boolean ggrd_grdtrack_interpolate_xy(double xin,double yin,
struct ggrd_gt *g,
double *value,
- unsigned char verbose)
+ ggrd_boolean verbose)
{
double x[3];
- unsigned char result;
+ ggrd_boolean result;
if(!g->init){ /* this will not necessarily work */
fprintf(stderr,"ggrd_grdtrack_interpolate_xy: error, g structure not initialized\n");
return FALSE;
@@ -393,14 +393,14 @@
struct GRD_HEADER **grd, /* pass as empty */
struct GMT_EDGEINFO **edgeinfo, /* pass as empty */
char *edgeinfo_string, /* -L type flags from GMT, can be empty */
- unsigned char *geographic_in, /* this is normally TRUE */
+ ggrd_boolean *geographic_in, /* this is normally TRUE */
int *pad, /* [4] array with padding (output) */
- unsigned char three_d, char *dfile, /* depth file name */
+ ggrd_boolean three_d, char *dfile, /* depth file name */
float **z, /* layers, pass as NULL */
int *nz, /* number of layers */
- unsigned char bilinear, /* linear/cubic? */
- unsigned char verbose,
- unsigned char change_depth_sign, /* change the
+ ggrd_boolean bilinear, /* linear/cubic? */
+ ggrd_boolean verbose,
+ ggrd_boolean change_depth_sign, /* change the
sign of the
depth
levels to go from depth (>0) to z (<0) */
@@ -409,10 +409,10 @@
int ggrd_grdtrack_init(double *west, double *east,double *south, double *north,
float **f,int *mm,char *grdfile,struct GRD_HEADER **grd,
struct GMT_EDGEINFO **edgeinfo,char *edgeinfo_string,
- unsigned char *geographic_in,int *pad,unsigned char three_d,
+ ggrd_boolean *geographic_in,int *pad,ggrd_boolean three_d,
char *dfile, float **z,int *nz,
- unsigned char bilinear,unsigned char verbose,
- unsigned char change_depth_sign)
+ ggrd_boolean bilinear,ggrd_boolean verbose,
+ ggrd_boolean change_depth_sign)
#endif
{
FILE *din;
@@ -699,8 +699,8 @@
*/
#ifdef USE_GMT4
-unsigned char ggrd_grdtrack_interpolate(double *in, /* lon/lat/z [2/3] in degrees/km */
- unsigned char three_d, /* use 3-D inetrpolation or 2-D? */
+ggrd_boolean ggrd_grdtrack_interpolate(double *in, /* lon/lat/z [2/3] in degrees/km */
+ ggrd_boolean three_d, /* use 3-D inetrpolation or 2-D? */
struct GRD_HEADER *grd, /* grd information */
float *f, /* data array */
struct GMT_EDGEINFO *edgeinfo, /* edge information */
@@ -708,11 +708,11 @@
float *z, /* depth layers */
int nz, /* number of depth layers */
double *value, /* output value */
- unsigned char verbose,
+ ggrd_boolean verbose,
struct GMT_BCR *bcr)
#else
-unsigned char ggrd_grdtrack_interpolate(double *in, /* lon/lat/z [2/3] in degrees/km */
- unsigned char three_d, /* use 3-D inetrpolation or 2-D? */
+ggrd_boolean ggrd_grdtrack_interpolate(double *in, /* lon/lat/z [2/3] in degrees/km */
+ ggrd_boolean three_d, /* use 3-D inetrpolation or 2-D? */
struct GRD_HEADER *grd, /* grd information */
float *f, /* data array */
struct GMT_EDGEINFO *edgeinfo, /* edge information */
@@ -720,12 +720,13 @@
float *z, /* depth layers */
int nz, /* number of depth layers */
double *value, /* output value */
- unsigned char verbose
+ ggrd_boolean verbose
)
#endif
{
int i1,i2;
double fac1,fac2,val1,val2;
+ static ggrd_boolean zwarned; /* this should move, leave for now */
/* If point is outside grd area,
shift it using periodicity or skip if not periodic. */
@@ -754,7 +755,7 @@
interpolate
*/
if(three_d){
- ggrd_gt_interpolate_z(in[2],z,nz,&i1,&i2,&fac1,&fac2,verbose);
+ ggrd_gt_interpolate_z(in[2],z,nz,&i1,&i2,&fac1,&fac2,verbose,&zwarned);
/*
we need these calls to reset the bcr.i and bcr.j counters
otherwise the interpolation routine would assume we have the same
@@ -820,10 +821,10 @@
returns error code
*/
-int ggrd_read_time_intervals(struct ggrd_t *thist,
- char *input_file,
- unsigned char read_thistory,
- unsigned char verbose)
+int ggrd_init_thist_from_file(struct ggrd_t *thist,
+ char *input_file,
+ ggrd_boolean read_thistory,
+ ggrd_boolean verbose)
{
FILE *in;
double ta,tb;
@@ -887,6 +888,7 @@
if(verbose)
fprintf(stderr,"ggrd_read_time_intervals: only one timestep\n");
}
+ thist->called = FALSE;
thist->init = TRUE;
return 0;
}
@@ -899,17 +901,17 @@
void ggrd_gt_interpolate_z(double z,float *za,int nz,
int *i1, int *i2,
double *fac1, double *fac2,
- unsigned char verbose)
+ ggrd_boolean verbose,
+ ggrd_boolean *zwarned)
{
int nzm1;
- static unsigned char warned=FALSE;
- if((!warned)&&verbose){
+ if((!(*zwarned))&&verbose){
if((z<za[0])||(z>za[nz-1])){
fprintf(stderr,"interpolate_z: WARNING: at least one z value extrapolated\n");
fprintf(stderr,"interpolate_z: zmin: %g z: %g zmax: %g\n",
za[0],z,za[nz-1]);
- warned = TRUE;
+ *zwarned = TRUE;
}
}
nzm1 = nz-1;
@@ -952,15 +954,7 @@
int *i1,int *i2,GGRD_CPREC *f1,GGRD_CPREC *f2,
GGRD_CPREC dxlimit)
{
- // dxlimit is the width of the transition between velocity stages
- static GGRD_CPREC xllimit,xrlimit;
- // these are the local, saved, quantities which are not determined anew
- // if the time doesn't change
- //
- static GGRD_CPREC f1_loc,f2_loc,time_old;
- static int ileft_loc,iright_loc,ntlim;
- static unsigned char called = FALSE;
GGRD_CPREC tloc,xll;
int i22,i;
if(!thist->init){
@@ -977,13 +971,13 @@
*f1 = 1.0;
*f2 = 0.0;
}else{ // more than one stage
- if(!called){
+ if(!thist->called){
/*
init loop
*/
/* check if dxlimit is larger than actual stage intervals */
- for(i=0;i<thist->nvtimes;i++){
+ for(i=0;i < thist->nvtimes;i++){
xll = thist->vtimes[i*3+2] - thist->vtimes[i*3];
if(dxlimit > xll){
fprintf(stderr,"inter_vel: adjusting transition width to %g from time intervals\n",
@@ -991,15 +985,15 @@
dxlimit = xll;
}
}
- ntlim = thist->nvtimes-1;
+ thist->ntlim = thist->nvtimes-1;
/* init some more factors */
- xllimit= -dxlimit/2.0;
- xrlimit= dxlimit/2.0;
- time_old = thist->vtimes[0] - 100; /* to make sure we get a
+ thist->xllimit= -dxlimit/2.0;
+ thist->xrlimit= dxlimit/2.0;
+ thist->time_old = thist->vtimes[0] - 100; /* to make sure we get a
first fix */
- called = TRUE;
+ thist->called = TRUE;
} /* end init */
- if(fabs(tloc - time_old)>5e-7){
+ if(fabs(tloc - thist->time_old)>5e-7){
/*
need to get new factors
*/
@@ -1029,19 +1023,19 @@
exit(-1);
}
}
- iright_loc=0; // right interval index
+ thist->iright_loc=0; // right interval index
i22=1; // right interval midpoint
// find the right interval such that its midpoint is larger or equal than time
while((tloc > thist->vtimes[i22]) &&
- (iright_loc < ntlim)){
- iright_loc++;
+ (thist->iright_loc < thist->ntlim)){
+ thist->iright_loc++;
i22 += 3;
}
- if(iright_loc == 0){
- iright_loc = 1;
+ if(thist->iright_loc == 0){
+ thist->iright_loc = 1;
i22 = 4;
}
- ileft_loc = iright_loc - 1; // left interval index
+ thist->ileft_loc = thist->iright_loc - 1; // left interval index
//
// distance from right boundary of left interval
// (=left boundary of right interval)
@@ -1055,29 +1049,29 @@
// vf1_loc and vf2_loc are the weights for velocities within the left and right
// intervals respectively
//
- if(xll < xllimit){ // xllimit should be 1-dx, dx~0.1
- f1_loc = 1.0;
- f2_loc = 0.0;
+ if(xll < thist->xllimit){ // xllimit should be 1-dx, dx~0.1
+ thist->f1_loc = 1.0;
+ thist->f2_loc = 0.0;
} else {
- if(xll > xrlimit){ // xrlimit should be 1+dx, dx~0.1
- f1_loc = 0.0;
- f2_loc = 1.0;
+ if(xll > thist->xrlimit){ // xrlimit should be 1+dx, dx~0.1
+ thist->f1_loc = 0.0;
+ thist->f2_loc = 1.0;
}else { // in between
- xll = (xll-xllimit)/dxlimit; // normalize by transition width
- f2_loc = ((1.0 - cos(xll * GGRD_PI))/2.0); // this goes from 0 to 1
+ xll = (xll-thist->xllimit)/dxlimit; // normalize by transition width
+ thist->f2_loc = ((1.0 - cos(xll * GGRD_PI))/2.0); // this goes from 0 to 1
// weight for left velocities
- f1_loc = 1.0 - f2_loc;
+ thist->f1_loc = 1.0 - thist->f2_loc;
}
}
//fprintf(stderr,"%g %g %g %i %i %g %g\n",time,f1_loc,f2_loc,ileft_loc,iright_loc,thist->vtimes[3*ileft_loc+2],thist->vtimes[3*iright_loc]);
- time_old = tloc;
+ thist->time_old = tloc;
}
// assign the local variables to the output variables
//
- *i1 = ileft_loc;
- *i2 = iright_loc;
- *f1 = f1_loc;
- *f2 = f2_loc;
+ *i1 = thist->ileft_loc;
+ *i2 = thist->iright_loc;
+ *f1 = thist->f1_loc;
+ *f2 = thist->f2_loc;
} // end ntimes>1 part
}
@@ -1085,26 +1079,25 @@
interpolate seafloor age at location vt,vp and time age
+this is a scalar interpolation that needs special care
+
*/
int interpolate_seafloor_ages(GGRD_CPREC xt, GGRD_CPREC xp,
GGRD_CPREC age,
struct ggrd_vel *v,
GGRD_CPREC *seafloor_age)
{
- static short int init = FALSE;
- static GGRD_CPREC old_age,old_f1,old_f2;
- static int old_left,old_right,ntlim;
int left, right,i;
GGRD_CPREC f1,f2;
double a1,a2;
- if(!init){
+ if(!v->sf_init){
/*
init the constants
*/
- ntlim = v->nage - 1;
+ v->sf_ntlim = v->nage - 1;
if(!v->nage){
fprintf(stderr,"interpolate_seafloor_ages: not initialized?!\n");
return -2;
@@ -1114,25 +1107,25 @@
fprintf(stderr,"interpolate_seafloor_ages: error: times need to be sorted monotnically increasing\n");
return -3;
}
- old_age = v->age_time[0] - 1000; /* so that we get interpolation
+ v->sf_old_age = v->age_time[0] - 1000; /* so that we get interpolation
factors */
- init = TRUE;
+ v->sf_init = TRUE;
}
/* check range */
- if((age < v->age_time[0]) || (age > v->age_time[ntlim])){
+ if((age < v->age_time[0]) || (age > v->age_time[v->sf_ntlim])){
fprintf(stderr,"interpolate_seafloor_ages: age: %g out of bounds [%g;%g]\n",
- age,v->age_time[0],v->age_time[ntlim]);
+ age,v->age_time[0],v->age_time[v->sf_ntlim]);
return -3;
}
- if(fabs(age-old_age) > 1e-8){
+ if(fabs(age- v->sf_old_age) > 1e-8){
/*
time interpolation
*/
right = 0;
- while((right < ntlim) && (v->age_time[right] < age))
+ while((right < v->sf_ntlim) && (v->age_time[right] < age))
right++;
if(right == 0)
right++;
@@ -1141,12 +1134,12 @@
f2 = (age - v->age_time[left])/(v->age_time[right]-v->age_time[left]);
f1 = 1.0-f2;
//fprintf(stderr,"sai: %g %g %g\t%g %g\n",v->age_time[left],age,v->age_time[right],f1,f2);
- old_age = age;
- old_left = left;old_right = right;
- old_f1 = f1;old_f2 = f2;
+ v->sf_old_age = age;
+ v->sf_old_left = left;v->sf_old_right = right;
+ v->sf_old_f1 = f1;v->sf_old_f2 = f2;
}else{ /* reuse time interpolation */
- left = old_left;right = old_right;
- f1 = old_f1;f2 = old_f2;
+ left = v->sf_old_left;right = v->sf_old_right;
+ f1 = v->sf_old_f1;f2 = v->sf_old_f2;
}
if(!ggrd_grdtrack_interpolate_tp((double)xt,(double)xp,
(v->ages+left),&a1,FALSE)){
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.h
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.h 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.h 2007-07-12 19:08:51 UTC (rev 7650)
@@ -18,28 +18,28 @@
wrappers
*/
-int ggrd_grdtrack_init_general(unsigned char ,char *, char *,char *,
- struct ggrd_gt *,unsigned char ,
- unsigned char);
-unsigned char ggrd_grdtrack_interpolate_rtp(double ,double ,double ,
+int ggrd_grdtrack_init_general(ggrd_boolean ,char *, char *,char *,
+ struct ggrd_gt *,ggrd_boolean ,
+ ggrd_boolean);
+ggrd_boolean ggrd_grdtrack_interpolate_rtp(double ,double ,double ,
struct ggrd_gt *,double *,
- unsigned char);
-unsigned char ggrd_grdtrack_interpolate_xyz(double ,double ,double ,
+ ggrd_boolean);
+ggrd_boolean ggrd_grdtrack_interpolate_xyz(double ,double ,double ,
struct ggrd_gt *,double *,
- unsigned char);
-unsigned char ggrd_grdtrack_interpolate_xy(double ,double ,
+ ggrd_boolean);
+ggrd_boolean ggrd_grdtrack_interpolate_xy(double ,double ,
struct ggrd_gt *,
double *,
- unsigned char );
-unsigned char ggrd_grdtrack_interpolate_tp(double ,double ,
+ ggrd_boolean );
+ggrd_boolean ggrd_grdtrack_interpolate_tp(double ,double ,
struct ggrd_gt *,
double *,
- unsigned char );
+ ggrd_boolean );
void ggrd_grdtrack_free_gstruc(struct ggrd_gt *);
-int ggrd_grdtrack_rescale(struct ggrd_gt *,unsigned char , unsigned char ,
- unsigned char ,double);
+int ggrd_grdtrack_rescale(struct ggrd_gt *,ggrd_boolean , ggrd_boolean ,
+ ggrd_boolean ,double);
/*
@@ -47,20 +47,20 @@
moderately external
*/
-int ggrd_read_time_intervals(struct ggrd_t *,char *,unsigned char ,unsigned char);
+int ggrd_init_thist_from_file(struct ggrd_t *,char *,ggrd_boolean ,ggrd_boolean);
int ggrd_read_vel_grids(struct ggrd_vel *, double, unsigned short, unsigned short, char *);
#ifdef USE_GMT4
/* GMT4.1.2 */
-unsigned char ggrd_grdtrack_interpolate(double *, unsigned char , struct GRD_HEADER *, float *,
- struct GMT_EDGEINFO *, int, float *, int , double *,unsigned char,
+ggrd_boolean ggrd_grdtrack_interpolate(double *, ggrd_boolean , struct GRD_HEADER *, float *,
+ struct GMT_EDGEINFO *, int, float *, int , double *,ggrd_boolean,
struct GMT_BCR *);
-int ggrd_grdtrack_init(double *, double *, double *, double *, float **, int *, char *, struct GRD_HEADER **, struct GMT_EDGEINFO **, char *, unsigned char *, int *, unsigned char, char *, float **, int *, unsigned char, unsigned char, unsigned char, struct GMT_BCR *);
+int ggrd_grdtrack_init(double *, double *, double *, double *, float **, int *, char *, struct GRD_HEADER **, struct GMT_EDGEINFO **, char *, ggrd_boolean *, int *, ggrd_boolean, char *, float **, int *, ggrd_boolean, ggrd_boolean, ggrd_boolean, struct GMT_BCR *);
#else
/* GMT 3.4.5 */
-unsigned char ggrd_grdtrack_interpolate(double *, unsigned char , struct GRD_HEADER *, float *,
- struct GMT_EDGEINFO *, int, float *, int , double *,unsigned char);
+ggrd_boolean ggrd_grdtrack_interpolate(double *, ggrd_boolean , struct GRD_HEADER *, float *,
+ struct GMT_EDGEINFO *, int, float *, int , double *,ggrd_boolean);
int ggrd_grdtrack_init(double *, double *,double *, double *, /* geographic bounds,
set all to zero to
@@ -72,13 +72,13 @@
char *, /* name, or prefix, of grd file with scalars */
struct GRD_HEADER **,
struct GMT_EDGEINFO **,
- char *,unsigned char *,
+ char *,ggrd_boolean *,
int *, /* [4] array with padding (output) */
- unsigned char _d, char *, /* depth file name */
+ ggrd_boolean _d, char *, /* depth file name */
float **, /* layers, pass as NULL */
int *, /* number of layers */
- unsigned char , /* linear/cubic? */
- unsigned char ,unsigned char);
+ ggrd_boolean , /* linear/cubic? */
+ ggrd_boolean ,ggrd_boolean);
void ggrd_gt_bcr_init_loc(void);
#endif
/*
@@ -89,7 +89,8 @@
void ggrd_gt_bcr_init_loc(void);
void ggrd_gt_interpolate_z(double,float *,int ,
- int *, int *, double *, double *,unsigned char); /* */
+ int *, int *, double *, double *,ggrd_boolean,
+ ggrd_boolean *); /* */
float ggrd_gt_rms(float *,int );
float ggrd_gt_mean(float *,int );
Modified: mc/1D/hc/trunk/ggrd_readgrds.c
===================================================================
--- mc/1D/hc/trunk/ggrd_readgrds.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/ggrd_readgrds.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -56,6 +56,21 @@
v->thist.init = FALSE;
v->velscale = 1.0;
v->rcmb = GGRD_RCMB_ND;
+ /* rlevels */
+ v->rlevels = NULL;
+ v->rl_warned = FALSE;
+ /* seafloor */
+ v->sf_init = FALSE;
+ /* thist structure */
+ v->thist.init = FALSE;
+ v->thist.called = FALSE;
+
+ /* interpolation structure */
+ v->vd.init = FALSE;
+ v->vd.reduce_r_stencil = FALSE;
+ v->vd.z_warned = FALSE;
+ v->vd.w_warned = FALSE;
+
}
/*
@@ -123,7 +138,7 @@
/*
if v->history is set, will look for different time intervals
*/
- ggrd_read_time_intervals(&v->thist,tfilename,v->history,verbose);
+ ggrd_init_thist_from_file(&v->thist,tfilename,v->history,verbose);
if(v->use_age){
/*
@@ -389,18 +404,18 @@
level);
ggrd_resort_and_check((v->vr+os1),fgrd,dgrd,v->n[HC_PHI],
v->n[HC_THETA],wraparound,1.0/v->velscale,
- v->read_gmt,TRUE,0.0);
+ v->read_gmt,TRUE,0.0,&v->vd.w_warned);
}else if((zero_boundary_vr)&&(v->rlevels[i] < v->rcmb)){
if(verbose)
fprintf(stderr,"ggrd_read_vel_grids: WARNING: assuming level %3i is at CMB and setting vr to zero\n",
level);
ggrd_resort_and_check((v->vr+os1),fgrd,dgrd,v->n[HC_PHI],
v->n[HC_THETA],wraparound,1.0/v->velscale,
- v->read_gmt,TRUE,0.0);
+ v->read_gmt,TRUE,0.0,&v->vd.w_warned);
}else
ggrd_resort_and_check((v->vr+os1),fgrd,dgrd,v->n[HC_PHI],
v->n[HC_THETA],wraparound,1.0/v->velscale,
- v->read_gmt,FALSE,ddummy);
+ v->read_gmt,FALSE,ddummy,&v->vd.w_warned);
hc_calc_mean_and_stddev((v->vr+os1),&ddummy,v->n[HC_TPPROD],
(mean+j),(std+j),(rms+j),FALSE,weighted,weights);
}else if(j == HC_THETA){
@@ -409,7 +424,7 @@
//
ggrd_resort_and_check((v->vt+os1),fgrd,dgrd,v->n[HC_PHI],
v->n[HC_THETA],wraparound,1.0/v->velscale,
- v->read_gmt,FALSE,ddummy);
+ v->read_gmt,FALSE,ddummy,&v->vd.w_warned);
hc_calc_mean_and_stddev((v->vt+os1),&ddummy,v->n[HC_TPPROD],
(mean+j),(std+j),(rms+j),FALSE,weighted,weights);
}else{
@@ -420,7 +435,7 @@
GGRD_PE("ggrd_read_vel_grds: index error");
ggrd_resort_and_check((v->vp+os1),fgrd,dgrd,v->n[HC_PHI],v->n[HC_THETA],
wraparound,1.0/v->velscale,
- v->read_gmt,FALSE,ddummy);
+ v->read_gmt,FALSE,ddummy,&v->vd.w_warned);
hc_calc_mean_and_stddev((v->vp+os1),&ddummy,v->n[HC_TPPROD],
(mean+j),(std+j),(rms+j),FALSE,weighted,weights);
//
@@ -476,19 +491,19 @@
int m, int n,hc_boolean wrap,
GGRD_CPREC factor,hc_boolean read_gmt,
hc_boolean set_to_constant,
- GGRD_CPREC constant)
+ GGRD_CPREC constant,
+ ggrd_boolean *warned)
{
int i,j,nm,os1,os2,boff;
- static hc_boolean warned = FALSE;
nm = m*n;
if(read_gmt){
// check for NaNs
for(i=0;i < nm;i++)
if(!finite(fb[i])){
fb[i] = 0.0;
- if(!warned){
+ if(!(*warned)){
fprintf(stderr,"WARNING: at least one NaN entry in the data has been replaced with zero\n");
- warned=TRUE;
+ *warned=TRUE;
}
}
}else{
@@ -496,9 +511,9 @@
for(i=0;i<nm;i++)
if(!finite(db[i])){
db[i]=0.0;
- if(!warned){
+ if(!(*warned)){
fprintf(stderr,"WARNING: at least one NaN entry in the data has been replaced with zero\n");
- warned=TRUE;
+ *warned=TRUE;
}
}
}
@@ -554,11 +569,11 @@
FILE *in;
int i;
GGRD_CPREC *rnew;
- hc_boolean warned = FALSE;
- in=hc_open(filename,"r","ggrd_read_depth_levels");
+
+ in = hc_open(filename,"r","ggrd_read_depth_levels");
/* set counters */
v->n[HC_R]=0;
- v->rlevels=(GGRD_CPREC *)malloc(sizeof(GGRD_CPREC));
+ v->rlevels=(GGRD_CPREC *)realloc(v->rlevels,sizeof(GGRD_CPREC));
if(!v->rlevels)
HC_MEMERROR("ggrd_read_depth_levels");
while(fscanf(in,HC_FLT_FORMAT,(v->rlevels + v->n[HC_R]))==1){
@@ -568,10 +583,10 @@
if(v->rlevels[v->n[HC_R]] < 0){
/* flip sign */
v->rlevels[v->n[HC_R]] = -v->rlevels[v->n[HC_R]];
- if((!warned) && (verbose)){
+ if((!v->rl_warned) && (verbose)){
fprintf(stderr,"ggrd_read_depth_levels: WARNING: flipping sign of depth levels in %s\n",
GGRD_DFILE);
- warned=TRUE;
+ v->rl_warned = TRUE;
}
}
/* radius of levels */
Modified: mc/1D/hc/trunk/ggrd_struc.h
===================================================================
--- mc/1D/hc/trunk/ggrd_struc.h 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/ggrd_struc.h 2007-07-12 19:08:51 UTC (rev 7650)
@@ -27,6 +27,12 @@
GGRD_CPREC tmin,tmax; /* range of times */
unsigned char init;
+
+ GGRD_CPREC xllimit,xrlimit;
+ GGRD_CPREC f1_loc,f2_loc,time_old;
+ int ileft_loc,iright_loc,ntlim;
+ unsigned char called;
+
};
/*
@@ -52,11 +58,21 @@
unsigned char init,
is_three; /* is it a 3-D set? */
+
+ double west,east,south,north;
+
#ifdef USE_GMT4
struct GMT_BCR bcr;
#endif
};
+/* velocity interpolation structure */
+struct ggrd_vip{
+ int ider[1+3*GGRD_MAX_IORDER],istencil[3],
+ ixtracer[3],old_order,orderp1,isshift[3];
+ unsigned char init,reduce_r_stencil,z_warned,w_warned;
+
+};
/*
@@ -78,6 +94,7 @@
history, /* time-dependent? */
use_age, /* use an additional age file */
read_gmt; /* read GMT grd files or binary format?*/
+ unsigned char rl_warned,vd_init,vd_reduce_r_stencil; /* */
int amode;
struct ggrd_gt *ages; /* for seafloor ages */
int nage; /* ntime for velo + 1 */
@@ -85,6 +102,11 @@
float age_bandlim; /* bandlim for age to decide on
continent
*/
+ struct ggrd_vip vd; /* velocity interpolation structure */
+ /* seafloor stuff */
+ unsigned short sf_init;
+ GGRD_CPREC sf_old_age,sf_old_f1,sf_old_f2;
+ int sf_old_left,sf_old_right,sf_ntlim;
};
#define GGRD_STRUC_INIT
#endif
Modified: mc/1D/hc/trunk/ggrd_velinterpol.c
===================================================================
--- mc/1D/hc/trunk/ggrd_velinterpol.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/ggrd_velinterpol.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -43,12 +43,9 @@
GGRD_CPREC w[GGRD_MAX_ORDERP1][GGRD_MAX_IORDERP1];
};
struct wgt weights[3];
- static int ider[1+3*GGRD_MAX_IORDER],istencil[3],
- ixtracer[3],old_order,orderp1,isshift[3];
+
- static hc_boolean init = FALSE, reduce_r_stencil = FALSE;
-
- if(!init){
+ if(!v->vd.init){
//
// do some checks if called for the first time
//
@@ -57,13 +54,13 @@
order,GGRD_MAX_ORDER);
return(-1);
}
- old_order = order;
- orderp1 = order+1;
+ v->vd.old_order = order;
+ v->vd.orderp1 = order+1;
if(v->n[HC_R] < order+1){
if(verbose)
fprintf(stderr,"ggrd_find_vel_and_der: WARNING: reducing r stencil to nl-1: %i\n",
v->n[HC_R] -1 );
- reduce_r_stencil = TRUE;
+ v->vd.reduce_r_stencil = TRUE;
}
if(v->n[HC_PHI] < order+1){
fprintf(stderr,"ggrd_find_vel_and_der: need at least four lon levels\n");
@@ -88,7 +85,7 @@
for(i=0;i < 3;i++){
if(i==HC_R)
- lorder = (reduce_r_stencil)?(v->n[HC_R]-1):(order);
+ lorder = (v->vd.reduce_r_stencil)?(v->n[HC_R]-1):(order);
else
lorder = order;
if(lorder < 0){
@@ -99,20 +96,20 @@
/* loop through dimensions */
/* all directions will be interpolated up to the
default order + 1 */
- istencil[i] = lorder + 1;
+ v->vd.istencil[i] = lorder + 1;
// these are offsets for the interpolation routine
- isshift[i] = (int)(istencil[i]/2.0);
+ v->vd.isshift[i] = (int)(v->vd.istencil[i]/2.0);
}
//
// this is for derivatives, initialize once as zeroes
//
for(i=0;i < 1+GGRD_MAX_IORDER*3;i++)
- ider[i] = 0.0;
- init = TRUE;
+ v->vd.ider[i] = 0.0;
+ v->vd.init = TRUE;
} /* end of init loop */
- if(order != old_order){
+ if(order != v->vd.old_order){
fprintf(stderr,"ggrd_find_vel_and_der: error: order (%i) shouldn't change, old: %i\n",
- order,old_order);
+ order,v->vd.old_order);
return(-1);
}
//
@@ -160,38 +157,38 @@
// RADIAL COMPONENT
//
// default order-1 interpolation for radial direction
- ixtracer[HC_R] = -1;
+ v->vd.ixtracer[HC_R] = -1;
ilim = v->n[HC_R] - 1;
i=0;
- while ((ixtracer[HC_R] == -1) && (i < ilim)){
+ while ((v->vd.ixtracer[HC_R] == -1) && (i < ilim)){
if(xloc[HC_R] <= v->rlevels[i]){
- ixtracer[HC_R] = i;
+ v->vd.ixtracer[HC_R] = i;
break;
}
i++;
}
- if(ixtracer[HC_R] == -1){ // no depth levels found, tracer is above surface
- ixtracer[HC_R] = ilim; // assign last layer, x_r should be corrected by the RK routines
+ if(v->vd.ixtracer[HC_R] == -1){ // no depth levels found, tracer is above surface
+ v->vd.ixtracer[HC_R] = ilim; // assign last layer, x_r should be corrected by the RK routines
}
//
// pick indices of grid points for the radial stencil
//
- for(i=0;i < istencil[HC_R];i++)
- igrid[HC_R][i] = ixtracer[HC_R] - isshift[HC_R] + i;
+ for(i=0;i < v->vd.istencil[HC_R];i++)
+ igrid[HC_R][i] = v->vd.ixtracer[HC_R] - v->vd.isshift[HC_R] + i;
//
// make sure all grid points exist
//
ishift = igrid[HC_R][0];
if(ishift < 0)
- for(i=0;i < istencil[HC_R];i++)
+ for(i=0;i < v->vd.istencil[HC_R];i++)
igrid[HC_R][i] -= ishift;
// same for upper limit
- ishift = igrid[HC_R][istencil[HC_R]-1] - ilim;
+ ishift = igrid[HC_R][v->vd.istencil[HC_R]-1] - ilim;
if(ishift > 0)
- for(i=0;i < istencil[HC_R];i++)
+ for(i=0;i < v->vd.istencil[HC_R];i++)
igrid[HC_R][i] -= ishift;
// find values of r for each grid point
- for(j=HC_R*orderp1,i=0;i < istencil[HC_R];i++)
+ for(j=HC_R*(v->vd.orderp1),i=0;i < v->vd.istencil[HC_R];i++)
grid[j+i] = v->rlevels[igrid[HC_R][i]];
@@ -199,37 +196,37 @@
// THETA COMPONENT
//
ilim = v->n[HC_THETA]-1;
- ixtracer[HC_THETA] = (int)(xloc[HC_THETA]/v->dtheta);
+ v->vd.ixtracer[HC_THETA] = (int)(xloc[HC_THETA]/v->dtheta);
// pick grid points
- for(i=0;i < istencil[HC_THETA];i++)
- igrid[HC_THETA][i] = ixtracer[HC_THETA] - isshift[HC_THETA]+i;
+ for(i=0;i < v->vd.istencil[HC_THETA];i++)
+ igrid[HC_THETA][i] = v->vd.ixtracer[HC_THETA] - v->vd.isshift[HC_THETA]+i;
//
// adust grid points to avoid wrap-around
//
ishift = igrid[HC_THETA][0];
if(ishift < 0){
- for(i=0;i < istencil[HC_THETA];i++)
+ for(i=0;i < v->vd.istencil[HC_THETA];i++)
igrid[HC_THETA][i] -= ishift;
}
// same for upper limit
- ishift = igrid[HC_THETA][istencil[HC_THETA]-1] - ilim;
+ ishift = igrid[HC_THETA][v->vd.istencil[HC_THETA]-1] - ilim;
if(ishift > 0)
- for(i=0;i < istencil[HC_THETA];i++)
+ for(i=0;i < v->vd.istencil[HC_THETA];i++)
igrid[HC_THETA][i] -= ishift;
//
// find values of theta: since given on dtheta/2 .... pi-dtheta/2
// theta_i = (i+0.5)*dtheta, i=0,1,...
//
- for(j=HC_THETA*orderp1,i=0;i < istencil[HC_THETA];i++)
+ for(j=HC_THETA*(v->vd.orderp1),i=0;i < v->vd.istencil[HC_THETA];i++)
grid[j+i] = (igrid[HC_THETA][i] + 0.5) * v->dtheta;
//
// now for phi
//
ilim = v->n[HC_PHI]-1;
- ixtracer[HC_PHI] = (int)(xloc[HC_PHI]/v->dphi+.5);
+ v->vd.ixtracer[HC_PHI] = (int)(xloc[HC_PHI]/v->dphi+.5);
//pick grid points
- for(i=0;i < istencil[HC_PHI];i++){
- igrid[HC_PHI][i] = ixtracer[HC_PHI] - isshift[HC_PHI] + i;
+ for(i=0;i < v->vd.istencil[HC_PHI];i++){
+ igrid[HC_PHI][i] = v->vd.ixtracer[HC_PHI] - v->vd.isshift[HC_PHI] + i;
//
// wrap around
//
@@ -241,7 +238,7 @@
//
// find values of phi. phi_i = i * dphi
//
- for(j=HC_PHI*orderp1,i=0;i < istencil[HC_PHI];i++)
+ for(j=HC_PHI*(v->vd.orderp1),i=0;i < v->vd.istencil[HC_PHI];i++)
grid[j+i] = igrid[HC_PHI][i] * v->dphi;
@@ -249,19 +246,19 @@
//
// check if all indices are ok
//
- for(i=0;i < istencil[HC_R];i++){
+ for(i=0;i < v->vd.istencil[HC_R];i++){
if((igrid[HC_R][i]< 0)||(igrid[HC_R][i] >= v->n[HC_R])){
fprintf(stderr,"ggrd_find_vel_and_der: row %i r index %i error\n",i,igrid[HC_R][i]);
return(-1);
}
}
- for(i=0;i < istencil[HC_THETA];i++){
+ for(i=0;i < v->vd.istencil[HC_THETA];i++){
if((igrid[HC_THETA][i] < 0) || (igrid[HC_THETA][i] >= v->n[HC_THETA])){
fprintf(stderr,"ggrd_find_vel_and_der: row %i theta index %i error \n",i,igrid[HC_THETA][i]);
return(-1);
}
}
- for(i=0;i < istencil[HC_PHI];i++){
+ for(i=0;i < v->vd.istencil[HC_PHI];i++){
if((igrid[HC_PHI][i] < 0)||(igrid[HC_PHI][i] >= v->n[HC_PHI])){
fprintf(stderr,"ggrd_find_vel_and_der: row %i phi index %i error\n",
i,igrid[HC_PHI][i]);
@@ -275,20 +272,20 @@
if(verbose >= 2){ /* debugging output */
fprintf(stderr,"ggrd_velinterpol: x={%g, %g, %g} [%i (%i), %i (%i), %i(%i)]\n",
- xloc[HC_R],xloc[HC_THETA],xloc[HC_PHI],ixtracer[HC_R],v->n[HC_R],
- ixtracer[HC_THETA],v->n[HC_THETA],
- ixtracer[HC_PHI],v->n[HC_PHI]);
+ xloc[HC_R],xloc[HC_THETA],xloc[HC_PHI],v->vd.ixtracer[HC_R],v->n[HC_R],
+ v->vd.ixtracer[HC_THETA],v->n[HC_THETA],
+ v->vd.ixtracer[HC_PHI],v->n[HC_PHI]);
for(i=0;i < 3;i++){
rnet = TWOPI;
fprintf(stderr,"ggrd_velinterpol: dim: %i:",i);
- for(j=0;j < istencil[i];j++){
- fprintf(stderr,"%.5f (%3i) ",grid[i*orderp1+j],igrid[i][j]);
+ for(j=0;j < v->vd.istencil[i];j++){
+ fprintf(stderr,"%.5f (%3i) ",grid[i*(v->vd.orderp1)+j],igrid[i][j]);
/* find min distance to stencil point */
- vrloc = fabs(grid[i*orderp1+j]-xloc[i]);
+ vrloc = fabs(grid[i*(v->vd.orderp1)+j]-xloc[i]);
if(vrloc < rnet){rnet=vrloc;k=j;}
}
fprintf(stderr,"\tms: %.4f(%i)\n",
- (GGRD_CPREC)k/(GGRD_CPREC)(istencil[i]-1),istencil[i]);
+ (GGRD_CPREC)k/(GGRD_CPREC)(v->vd.istencil[i]-1),v->vd.istencil[i]);
}
}
#endif
@@ -298,7 +295,7 @@
// compute all the weights for each stencil
//
for(i=0;i < 3;i++){ /* loop through spatial dimension */
- ggrd_weights(xloc[i],(grid+i*orderp1),istencil[i],iorder,weights[i].w);
+ ggrd_weights(xloc[i],(grid+i*(v->vd.orderp1)),v->vd.istencil[i],iorder,weights[i].w);
}
//
// first calculate velocities only (idindex=0) or vel and derivatives
@@ -309,9 +306,9 @@
first velocities
*/
dvr[0] = dvtheta[0] = dvphi[0] = 0.0;
- for(i=0;i < istencil[HC_R];i++){ // radial
- for(j=0; j < istencil[HC_THETA];j++){ // theta
- for(k=0; k < istencil[HC_PHI];k++){ // phi
+ for(i=0;i < v->vd.istencil[HC_R];i++){ // radial
+ for(j=0; j < v->vd.istencil[HC_THETA];j++){ // theta
+ for(k=0; k < v->vd.istencil[HC_PHI];k++){ // phi
rnet = weights[HC_R].w[i][0];
rnet *= weights[HC_THETA].w[j][0];
rnet *= weights[HC_PHI].w[k][0];
@@ -336,13 +333,13 @@
// this is the derivative yes/no array
// set once, and switch off again below//
//
- ider[m] = 1;
- for(i=0;i < istencil[HC_R];i++){
- for(j=0; j < istencil[HC_THETA];j++){
- for(k=0; k < istencil[HC_PHI];k++){
- rnet = weights[HC_R].w[i][ider[HC_R+1]];
- rnet *= weights[HC_THETA].w[j][ider[HC_THETA+1]];
- rnet *= weights[HC_PHI].w[k][ider[HC_PHI+1]];
+ v->vd.ider[m] = 1;
+ for(i=0;i < v->vd.istencil[HC_R];i++){
+ for(j=0; j < v->vd.istencil[HC_THETA];j++){
+ for(k=0; k < v->vd.istencil[HC_PHI];k++){
+ rnet = weights[HC_R].w[i][v->vd.ider[HC_R+1]];
+ rnet *= weights[HC_THETA].w[j][v->vd.ider[HC_THETA+1]];
+ rnet *= weights[HC_PHI].w[k][v->vd.ider[HC_PHI+1]];
index = igrid[HC_R][i] * v->n[HC_TPPROD] +
igrid[HC_THETA][j] * v->n[HC_PHI] +
igrid[HC_PHI][k];
@@ -354,7 +351,7 @@
}
}
}
- ider[m]=0; // reset derivative switxh to zero
+ v->vd.ider[m]=0; // reset derivative switxh to zero
}
/* succesful return */
return 0;
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hc.h 2007-07-12 19:08:51 UTC (rev 7650)
@@ -92,6 +92,22 @@
HC_PREC u[6][4];
};
/*
+ for polsol and solve
+*/
+struct hc_ps{
+ /* scaling factors will only be computed once */
+ int ncalled;
+ double rho_scale;
+ double alpha, beta, geoid_factor;
+ hc_boolean rho_init,
+ prop_params_init, /* parameters for propagator computation */
+ abg_init , /* alpha, beta factors */
+ prop_mats_init; /* will be true only if save_prop_mats is
+ requested */
+ /* for solve */
+ hc_boolean tor_init, pol_init;
+};
+/*
parameter structure to allow for settings that are specific to the
@@ -186,7 +202,7 @@
*/
-
+ struct hc_ps psp;
hc_boolean save_solution; /* memory intensive speedup in poloidal
solution by saving propagator matrices
this will also keep the toroidal solution
@@ -212,7 +228,7 @@
struct sh_lms *tor_sol; /*
toroidal solution
*/
- hc_boolean initialized; /* logic flag */
+ hc_boolean initialized,const_init,visc_init,dens_init,pvel_init; /* logic flags */
/* sqrt(l(l+1)) and 1/lfac factors */
HC_PREC *lfac,*ilfac;
int lfac_init;
@@ -237,7 +253,7 @@
HC_PREC stress_scale; /* to go to MPa */
HC_PREC r_cmb; /* radius of CMB */
/* Legendre functions */
- double *plm;
+ SH_RICK_PREC *plm;
/* more logic flags */
my_boolean spectral_solution_computed, spatial_solution_computed;
Modified: mc/1D/hc/trunk/hc2gmt
===================================================================
--- mc/1D/hc/trunk/hc2gmt 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hc2gmt 2007-07-12 19:08:51 UTC (rev 7650)
@@ -1,6 +1,6 @@
#!/bin/bash
#
-# convert HC output in binary format to GMT grd files
+# convert HC velocity output in binary format (as produced by -px) to GMT grd files for vr, vt, vp
#
dfile=vdepth.dat # depth level file
sfile=vsol # prefix of solution files
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2007-07-12 19:08:51 UTC (rev 7650)
@@ -8,8 +8,8 @@
void ggrd_grdtrack_free_gstruc(struct ggrd_gt *);
void ggrd_find_spherical_vel_from_rigid_cart_rot(double *, double *, double *, double *, double *);
void ggrd_print_layer_avg(float *, float *, int, int, FILE *);
-int ggrd_read_time_intervals(struct ggrd_t *, char *, unsigned char, unsigned char);
-void ggrd_gt_interpolate_z(double, float *, int, int *, int *, double *, double *, unsigned char);
+int ggrd_init_thist_from_file(struct ggrd_t *, char *, unsigned char, unsigned char);
+void ggrd_gt_interpolate_z(double, float *, int, int *, int *, double *, double *, unsigned char, unsigned char *);
void ggrd_interpol_time(double, struct ggrd_t *, int *, int *, double *, double *, double);
int interpolate_seafloor_ages(double, double, double, struct ggrd_vel *, double *);
float ggrd_gt_rms(float *, int);
@@ -17,7 +17,7 @@
/* ggrd_readgrds.c */
void ggrd_init_vstruc(struct ggrd_vel *);
int ggrd_read_vel_grids(struct ggrd_vel *, double, unsigned short, unsigned short, char *);
-void ggrd_resort_and_check(double *, float *, double *, int, int, unsigned short, double, unsigned short, unsigned short, double);
+void ggrd_resort_and_check(double *, float *, double *, int, int, unsigned short, double, unsigned short, unsigned short, double, unsigned char *);
void ggrd_read_depth_levels(struct ggrd_vel *, int **, char *, unsigned short);
/* ggrd_test.c */
/* ggrd_velinterpol.c */
@@ -28,6 +28,7 @@
/* hc_init.c */
void hc_init_parameters(struct hc_parameters *);
void hc_struc_init(struct hcs **);
+void hc_init_polsol_struct(struct hc_ps *);
void hc_init_main(struct hcs *, int, struct hc_parameters *);
void hc_init_constants(struct hcs *, double, char *, unsigned short);
void hc_handle_command_line(int, char **, struct hc_parameters *);
@@ -37,6 +38,7 @@
void hc_assign_plate_velocities(struct hcs *, int, char *, unsigned short, int, unsigned short, unsigned short);
void hc_init_l_factors(struct hcs *, int);
void hc_get_blank_expansions(struct sh_lms **, int, int, char *);
+void hc_struc_free(struct hcs **);
/* hc_input.c */
void hc_read_sh_solution(struct hcs *, struct sh_lms **, FILE *, unsigned short, unsigned short);
/* hc_matrix.c */
@@ -108,15 +110,18 @@
void rick_realft_nr(float *, int, int);
void rick_four1_nr(float *, int, int);
/* rick_sh_c.c */
-void rick_compute_allplm(int, int, double *, double *, struct rick_module *);
-void rick_pix2ang(int, int, double *, double *, struct rick_module *);
+void rick_compute_allplm(int, int, float *, float *, struct rick_module *);
+void rick_compute_allplm_irreg(int, int, float *, float *, struct rick_module *, float *, int);
+void rick_pix2ang(int, int, float *, float *, struct rick_module *);
void rick_shc2d(float *, float *, int, int, float *, float *, struct rick_module *);
-void rick_shc2d_pre(float *, float *, int, double *, double *, int, float *, float *, struct rick_module *);
+void rick_shc2d_irreg(float *, float *, int, int, float *, float *, struct rick_module *, float *, int, float *, int, unsigned short);
+void rick_shc2d_pre(float *, float *, int, float *, float *, int, float *, float *, struct rick_module *);
+void rick_shc2d_pre_irreg(float *, float *, int, float *, float *, int, float *, float *, struct rick_module *, float *, int, float *, int, unsigned short);
void rick_shd2c(float *, float *, int, int, float *, float *, struct rick_module *);
-void rick_shd2c_pre(float *, float *, int, double *, double *, int, float *, float *, struct rick_module *);
-void rick_init(int, int, int *, int *, int *, struct rick_module *);
+void rick_shd2c_pre(float *, float *, int, float *, float *, int, float *, float *, struct rick_module *);
+void rick_init(int, int, int *, int *, int *, struct rick_module *, unsigned short);
void rick_free_module(struct rick_module *, int);
-void rick_plmbar1(double *, double *, int, int, float, struct rick_module *);
+void rick_plmbar1(float *, float *, int, int, float, struct rick_module *);
void rick_gauleg(float, float, float *, float *, int);
/* sh_ana.c */
/* shana_sh.c */
@@ -130,8 +135,8 @@
void shana_free_module(struct shana_module *, int);
void shana_plmbar1(double *, double *, int, int, double, struct shana_module *);
/* sh_exp.c */
-void sh_allocate_and_init(struct sh_lms **, int, int, int, int, unsigned short);
-void sh_init_expansion(struct sh_lms *, int, int, int, unsigned short);
+void sh_allocate_and_init(struct sh_lms **, int, int, int, int, unsigned short, unsigned short);
+void sh_init_expansion(struct sh_lms *, int, int, int, unsigned short, unsigned short);
void sh_free_expansion(struct sh_lms *, int);
void sh_clear_alm(struct sh_lms *);
double sh_total_power(struct sh_lms *);
@@ -143,12 +148,15 @@
void sh_print_nonzero_coeff(struct sh_lms *, FILE *);
void sh_read_spatial_data_from_file(struct sh_lms *, FILE *, unsigned short, int, float *, float *);
void sh_compute_spatial_basis(struct sh_lms *, FILE *, unsigned short, float, float **, int, unsigned short);
-void sh_compute_spectral(float *, int, unsigned short, double **, struct sh_lms *, unsigned short);
-void sh_compute_spatial(struct sh_lms *, int, unsigned short, double **, float *, unsigned short);
+void sh_compute_spectral(float *, int, unsigned short, float **, struct sh_lms *, unsigned short);
+void sh_compute_spatial(struct sh_lms *, int, unsigned short, float **, float *, unsigned short);
+void sh_compute_spatial_irreg(struct sh_lms *, int, unsigned short, float **, float *, int, float *, int, float *, unsigned short, unsigned short);
void sh_exp_type_error(char *, struct sh_lms *);
-void sh_print_plm(double *, int, int, int, FILE *);
+void sh_print_plm(float *, int, int, int, FILE *);
void sh_print_spatial_data_to_file(struct sh_lms *, int, float *, unsigned short, float, FILE *);
-void sh_compute_plm(struct sh_lms *, int, double **, unsigned short);
+void sh_print_irreg_spatial_data_to_file(struct sh_lms *, int, float *, unsigned short, float, float *, int, float *, int, FILE *);
+void sh_compute_plm(struct sh_lms *, int, float **, unsigned short);
+void sh_compute_plm_irreg(struct sh_lms *, int, float **, unsigned short, float *, int);
void sh_get_coeff(struct sh_lms *, int, int, int, unsigned short, double *);
void sh_write_coeff(struct sh_lms *, int, int, int, unsigned short, double *);
void sh_aexp_equals_bexp_coeff(struct sh_lms *, struct sh_lms *);
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hc_init.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -52,10 +52,14 @@
*/
void hc_struc_init(struct hcs **hc)
{
-
+ /* this will take care of all flags and such */
*hc = (struct hcs *)calloc(1,sizeof(struct hcs ));
if(!(*hc))
HC_MEMERROR("hc_struc_init: hc");
+ /* just to be sure */
+ (*hc)->initialized = (*hc)->const_init = (*hc)->visc_init =
+ (*hc)->dens_init = (*hc)->pvel_init = FALSE;
+ hc_init_polsol_struct(&((*hc)->psp));
/*
assign NULL pointers to allow reallocating
*/
@@ -66,6 +70,18 @@
(*hc)->plm = NULL;
(*hc)->prem_init = FALSE;
}
+
+void hc_init_polsol_struct(struct hc_ps *psp)
+{
+ psp->ncalled = 0;
+ /* scaling factors will only be computed once */
+ psp->rho_scale = 1.0;
+ psp->rho_init = FALSE;
+ psp->prop_params_init = FALSE; /* parameters for propagator computation */
+ psp->abg_init = FALSE; /* alpha, beta factors */
+ psp->prop_mats_init = FALSE; /* will be true only if save_prop_mats is */
+ psp->tor_init = psp->pol_init = FALSE;
+}
/*
initialize all variables, after initializing the parameters
@@ -152,9 +168,8 @@
void hc_init_constants(struct hcs *hc, HC_PREC dens_anom_scale,
char *prem_filename,hc_boolean verbose)
{
- static hc_boolean init=FALSE;
int ec;
- if(init)
+ if(hc->const_init)
HC_ERROR("hc_init_constants","why would you call this routine twice?")
if(!hc->prem_init){
/* PREM constants */
@@ -220,7 +235,7 @@
(hc->timesc * 1e6);
- init = TRUE;
+ hc->const_init = TRUE;
}
/*
@@ -348,7 +363,6 @@
FILE *in;
char fstring[100];
HC_PREC mean,mweight,rold,mws;
- static hc_boolean init=FALSE;
switch(mode){
case HC_INIT_FROM_FILE:
/*
@@ -356,7 +370,7 @@
init from file part
*/
- if(init)
+ if(hc->visc_init)
HC_ERROR("hc_assign_viscosity","viscosity already read from file, really read again?");
/*
read viscosity structure from file
@@ -430,7 +444,7 @@
HC_ERROR("hc_assign_viscosity","mode undefined");
break;
}
- init = TRUE;
+ hc->visc_init = TRUE;
}
/*
@@ -477,18 +491,17 @@
to 0.01 for percent input,
for instance
*/
- static hc_boolean init=FALSE;
hc->compressible = compressible;
- if(init) /* clear old expansions, if
- already initialized */
+ if(hc->dens_init) /* clear old expansions, if
+ already initialized */
sh_free_expansion(hc->dens_anom,hc->inho);
/* get PREM model, if not initialized */
if(!hc->prem_init)
HC_ERROR("hc_assign_density","assign 1-D reference model (PREM) first");
switch(mode){
case HC_INIT_FROM_FILE:
- if(init)
+ if(hc->dens_init)
HC_ERROR("hc_assign_density","really read dens anomalies again from file?");
/*
@@ -573,11 +586,11 @@
hc_get_blank_expansions(&hc->dens_anom,hc->inho+1,hc->inho,
"hc_assign_density");
/*
- initialize expansion
+ initialize expansion on irregular grid
*/
sh_init_expansion((hc->dens_anom+hc->inho),
(nominal_lmax > lmax) ? (nominal_lmax):(lmax),
- hc->sh_type,1,verbose);
+ hc->sh_type,1,verbose,FALSE);
/*
read parameters and scale (put possible depth dependence of
@@ -599,7 +612,7 @@
HC_ERROR("hc_assign_density","mode undefined");
break;
}
- if((!init)||(layer_structure_changed)){
+ if((!hc->dens_init)||(layer_structure_changed)){
/*
assign the other radii, nrad + 2
@@ -649,7 +662,7 @@
sqrt(sh_total_power((hc->dens_anom+i))));
free(dbot);free(dtop);
} /* end layer structure part */
- init = TRUE;
+ hc->dens_init = TRUE;
}
/*
@@ -694,13 +707,12 @@
hc_boolean pvel_in_binary,
hc_boolean verbose)
{
- static hc_boolean init = FALSE;
int type,shps,ilayer,nset,ivec;
HC_PREC zlabel,vfac[2],t10[2],t11[2];
FILE *in;
/* scale to go from cm/yr to internal scale */
vfac[0] = vfac[1] = 1.0/hc->vel_scale;
- if(init)
+ if(hc->pvel_init)
HC_ERROR("hc_assign_plate_velocities","what to do if called twice?");
if(!vel_bc_zero){
/*
@@ -737,10 +749,10 @@
exit(-1);
}
/*
- initialize expansion
+ initialize expansion using irregular grid
*/
- sh_init_expansion((hc->pvel+0),lmax,hc->sh_type,1,verbose);
- sh_init_expansion((hc->pvel+1),lmax,hc->sh_type,1,verbose);
+ sh_init_expansion((hc->pvel+0),lmax,hc->sh_type,1,verbose, FALSE);
+ sh_init_expansion((hc->pvel+1),lmax,hc->sh_type,1,verbose, FALSE);
/*
read in expansions, convert to internal format from
physical
@@ -774,20 +786,21 @@
/*
initialize with zeroes
*/
- if(init){
+ if(hc->pvel_init){
sh_clear_alm(hc->pvel);
sh_clear_alm((hc->pvel+1));
}else{
+ /* use irregular grid */
sh_init_expansion(hc->pvel,lmax,hc->sh_type,
- 1,verbose);
+ 1,verbose,FALSE);
sh_init_expansion((hc->pvel+1),lmax,hc->sh_type,
- 1,verbose);
+ 1,verbose,FALSE);
}
if(verbose)
fprintf(stderr,"hc_assign_plate_velocities: using no-slip surface BC, lmax %i\n",
lmax);
}
- init = TRUE;
+ hc->pvel_init = TRUE;
}
/*
@@ -859,3 +872,23 @@
memcpy((*expansion+nold),tmpzero,ngrow*sizeof(struct sh_lms));
free(tmpzero);
}
+
+/*
+
+
+be more careful with freeing
+
+
+ */
+void hc_struc_free(struct hcs **hc)
+{
+ free((*hc)->visc);
+ free((*hc)->rvisc);
+ free((*hc)->visc);
+
+ free((*hc)->qwrite);
+
+ sh_free_expansion((*hc)->dens_anom,1);
+
+ free(*hc);
+}
Modified: mc/1D/hc/trunk/hc_input.c
===================================================================
--- mc/1D/hc/trunk/hc_input.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hc_input.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -36,8 +36,8 @@
if(!(*sol))
HC_MEMERROR("hc_read_sh_solution: sol");
hc->sh_type = type;
- for(i=0;i < nsol;i++) /* init expansions */
- sh_init_expansion((*sol+i),lmax,hc->sh_type,1,verbose);
+ for(i=0;i < nsol;i++) /* init expansions on irregular grid */
+ sh_init_expansion((*sol+i),lmax,hc->sh_type,1,verbose,FALSE);
hc->nrad = nset - 2;
hc_vecrealloc(&hc->r,nset,"hc_read_sh_solution");
}
Modified: mc/1D/hc/trunk/hc_misc.c
===================================================================
--- mc/1D/hc/trunk/hc_misc.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hc_misc.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -156,7 +156,7 @@
void hc_get_flt_frmt_string(char *string, int n,
hc_boolean append)
{
- static hc_boolean init=FALSE;
+ static hc_boolean init=FALSE; /* that's OK, multiple instances calling are fine */
static char type[2];
int i;
if(!init){
Modified: mc/1D/hc/trunk/hc_polsol.c
===================================================================
--- mc/1D/hc/trunk/hc_polsol.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hc_polsol.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -197,20 +197,10 @@
HC_PREC *xprem;
double *b,du1,du2,el,rnext,drho,dadd;
double amat[3][3],bvec[3],u[4],poten[2],unew[4],potnew[2],clm[2];
- static int ncalled = 0;
- /* scaling factors will only be computed once */
- static double rho_scale = 1.0;
- static double alpha, beta, geoid_factor;
/*
structures which hold u[6][4] type arrays
*/
struct hc_sm cmb, *u3;
- /* init flags */
- static hc_boolean rho_init = FALSE,
- prop_params_init = FALSE, /* parameters for propagator computation */
- abg_init = FALSE, /* alpha, beta factors */
- prop_mats_init = FALSE; /* will be true only if save_prop_mats is
- requested */
hc_boolean qvis,qinho,hit;
/*
define a few offset and size pointers
@@ -233,7 +223,7 @@
/*
check if still same general number of layers
*/
- if((prop_params_init)&&((nradp2 != hc->nradp2)||
+ if((hc->psp.prop_params_init)&&((nradp2 != hc->nradp2)||
(inho2 != hc->inho2)||
(nvisp1 != hc->nvisp1))){
HC_ERROR("hc_polsol","layer structure changed from last call");
@@ -252,7 +242,7 @@
/*
propagators saved
*/
- if(!prop_mats_init){
+ if(!hc->psp.prop_mats_init){
/*
we will be saving all propagator matrices. this makes sense if
@@ -273,26 +263,26 @@
hc_dvecalloc(&hc->props,prop_s1,"hc_polsol");
hc_dvecalloc(&hc->ppots,prop_s2,"hc_polsol");
}
- if(!abg_init){
+ if(!hc->psp.abg_init){
//
// SET alpha, beta and geoid factors
//
- alpha = rho_scale * (hc->re*10.) * hc->gacc / hc->visnor; /* */
- alpha *= ONEEIGHTYOVERPI * hc->secyr * hc->timesc; /* */
+ hc->psp.alpha = hc->psp.rho_scale * (hc->re*10.) * hc->gacc / hc->visnor; /* */
+ hc->psp.alpha *= ONEEIGHTYOVERPI * hc->secyr * hc->timesc; /* */
//
- beta = -4.0 * HC_PI * (hc->g*1e3) * (hc->re*1e2) / hc->gacc;
+ hc->psp.beta = -4.0 * HC_PI * (hc->g*1e3) * (hc->re*1e2) / hc->gacc;
if(verbose)
fprintf(stderr,"hc_polsol: alpha: %.8f beta: %.8f\n",
- alpha,beta);
+ hc->psp.alpha,hc->psp.beta);
/*
geoid scaling factor hc->gacc shoud be grav[nprops] for
compressibility
*/
- geoid_factor = HC_PI * 10.0* hc->visnor/180./hc->secyr/hc->gacc/1.e8;
- abg_init = TRUE;
+ hc->psp.geoid_factor = HC_PI * 10.0* hc->visnor/180./hc->secyr/hc->gacc/1.e8;
+ hc->psp.abg_init = TRUE;
}
- if((!prop_params_init) || (viscosity_or_layer_changed)){
+ if((!hc->psp.prop_params_init) || (viscosity_or_layer_changed)){
/*
intialize arrays that depend on viscosity and density layer spacing
@@ -308,7 +298,7 @@
// 5) IF AT RDEN(NIH) EVALUATE DEN, INCREMENT NIH
// 6) IF AT RAD(I) QWRITE = TRUE, INCREMENT I
//
- if(!prop_params_init){
+ if(!hc->psp.prop_params_init){
if(verbose)
fprintf(stderr,"hc_polsol: initializing for %i v layers and %i dens layers\n",
nrad,inho);
@@ -402,7 +392,7 @@
//
// IF RDEN, EVALUATE DEN, INCREMENT NIH
//
- hc->den[hc->nprops-1] = dfact[nih] * hc->rden[nih] * hc->rden[nih] * alpha;
+ hc->den[hc->nprops-1] = dfact[nih] * hc->rden[nih] * hc->rden[nih] * hc->psp.alpha;
nih++;
}
}
@@ -424,7 +414,7 @@
*/
if(verbose >= 3){
- if(prop_params_init)
+ if(hc->psp.prop_params_init)
fprintf(stderr,"hc_polsol: using old parameters: %i v layers and %i dens layers\n",
nrad,inho);
for(i=i2=0;i < hc->nprops+1;i++){
@@ -435,7 +425,7 @@
hc->pvisc[i],hc->den[i],i2,inho);
}
}
- if(!rho_init){
+ if(!hc->psp.rho_init){
/*
initialize the density factors, for incompressible, those
are all constant, else from PREM
@@ -478,9 +468,9 @@
hc->rho_zero[i] = hc->avg_den_mantle;
}
hc->rho_zero[hc->nprops+1] = 0.0;
- rho_init = TRUE;
+ hc->psp.rho_init = TRUE;
} /* end rho init */
- prop_params_init = TRUE;
+ hc->psp.prop_params_init = TRUE;
/*
end of the propagator factor section, this will only get executed
@@ -499,7 +489,7 @@
//
if(verbose)
fprintf(stderr,"hc_polsol: ncalled: %5i for lmax: %i dens lmax: %i, visc or layer %s changed\n",
- ncalled,pol_sol[0].lmax,dens_anom->lmax,
+ hc->psp.ncalled,pol_sol[0].lmax,dens_anom->lmax,
((viscosity_or_layer_changed)?(""):("not")));
if(free_slip) /* select which components of pol solvec to
use */
@@ -518,7 +508,7 @@
*/
el = (HC_PREC)l;
- if((!save_prop_mats) || (!prop_mats_init)){
+ if((!save_prop_mats) || (!hc->psp.prop_mats_init)){
//
// get all propagators now, as they only depend on l
//
@@ -610,7 +600,7 @@
if we were allowing for compressibility, would multi with
hc->grav[i]/hc->grav here (beta incorporates 1/grav0)
*/
- poten[1] += beta * hc->rprops[0] *
+ poten[1] += hc->psp.beta * hc->rprops[0] *
(u[2] - (hc->rho_zero[0] - hc->rho_zero[-1]) *
poten[0]);
@@ -665,7 +655,7 @@
//
// fprintf(stderr,"%15.5e %15.5e %15.5e %15.5e\n",
// beta, hc->den[i], b[ninho],hc->rden[ninho]);
- poten[1] += beta * dadd * hc->rden[ninho];
+ poten[1] += hc->psp.beta * dadd * hc->rden[ninho];
}
//
// Changes due to radial density variations
@@ -685,7 +675,7 @@
if ((jpb < npb)&&(hc->rprops[ip1] > rpb[jpb] - 0.0001)){
if (ibv == 0) {
u[2] -= fpb[jpb] * b[ninho];
- poten[1] -= beta * rpb[jpb] * fpb[jpb] * b[ninho] * rho_scale;
+ poten[1] -= hc->psp.beta * rpb[jpb] * fpb[jpb] * b[ninho] * hc->psp.rho_scale;
}
jpb++;
}
@@ -717,7 +707,7 @@
// is proportional to total surface elevation (not minus equipotential
// surface)
//
- poten[1] -= beta * hc->rprops[hc->nprops] *
+ poten[1] -= hc->psp.beta * hc->rprops[hc->nprops] *
(u[2] - hc->rho_zero[hc->nprops] * poten[0]);
nl = ilayer;
u3[nl].u[5][ibv] = poten[1];
@@ -801,7 +791,7 @@
} /* end l loop */
if(save_prop_mats)
/* only now can we set the propagator matrix storage scheme to TRUE */
- prop_mats_init = TRUE;
+ hc->psp.prop_mats_init = TRUE;
if(verbose)
fprintf(stderr,"hc_polsol: assigned nl: %i nprop: %i nrad: %i layers\n",
nl,hc->nprops,nrad);
@@ -831,12 +821,12 @@
for(m=0;m <= l;m++){
if (m != 0){
sh_get_coeff((pol_sol+os),l,m,2,FALSE,clm); /* internal convention */
- clm[0] *= geoid_factor;
- clm[1] *= geoid_factor;
+ clm[0] *= hc->psp.geoid_factor;
+ clm[1] *= hc->psp.geoid_factor;
sh_write_coeff(geoid,l,m,2,FALSE,clm);
}else{ /* m == 0 */
sh_get_coeff((pol_sol+os),l,m,0,FALSE,clm);
- clm[0] *= geoid_factor;
+ clm[0] *= hc->psp.geoid_factor;
sh_write_coeff(geoid,l,m,0,FALSE,clm);
}
}
@@ -856,7 +846,7 @@
free(hc->props);free(hc->ppots);
}
/* all others should be saved */
- ncalled++;
+ hc->psp.ncalled++;
if(verbose)
fprintf(stderr,"hc_polsol: done\n");
}
Modified: mc/1D/hc/trunk/hc_solve.c
===================================================================
--- mc/1D/hc/trunk/hc_solve.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hc_solve.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -52,9 +52,7 @@
compare with Benhard's
densub densol output */
double *tvec;
- static hc_boolean
- tor_init = FALSE,
- pol_init = FALSE;
+
if(!hc->initialized)
HC_ERROR("hc_solve","hc structure not initialized");
if((!free_slip) && (hc->pvel[0].lmax < hc->dens_anom[0].lmax)){
@@ -74,13 +72,13 @@
initialize a bunch of expansions for the poloidal solution
*/
nsh_pol = 6 * (hc->nrad+2); /* u[4] plus poten[2] */
- if((!pol_init)||(!hc->save_solution)){
+ if((!hc->psp.pol_init)||(!hc->save_solution)){
/* room for pol solution */
sh_allocate_and_init(&hc->pol_sol,nsh_pol,
hc->dens_anom[0].lmax,hc->sh_type,
- 0,verbose);
+ 0,verbose,FALSE); /* irregular grid */
}
- if((!hc->save_solution) || (!pol_init) || viscosity_or_layer_changed ||
+ if((!hc->save_solution) || (!hc->psp.pol_init) || viscosity_or_layer_changed ||
dens_anom_changed || ((!free_slip) && (plate_vel_changed))){
/*
@@ -110,12 +108,12 @@
solve toroidal part only for no-slip surface boundary condition
*/
- if((!tor_init)||(!hc->save_solution)){
+ if((!hc->psp.tor_init)||(!hc->save_solution)){
nsh_tor = 2 * (hc->nrad+2);
sh_allocate_and_init(&hc->tor_sol,nsh_tor,hc->pvel[1].lmax,
- hc->sh_type,0,verbose);
+ hc->sh_type,0,verbose,FALSE); /* irregular grid */
}
- if((!tor_init) || viscosity_or_layer_changed || plate_vel_changed ||
+ if((!hc->psp.tor_init) || viscosity_or_layer_changed || plate_vel_changed ||
(!hc->save_solution)){
/*
if we are not saving solutions, or the velocities or viscosities
@@ -176,8 +174,8 @@
if(!free_slip)
sh_free_expansion(hc->tor_sol,nsh_tor);
}
- pol_init = TRUE;
- tor_init = TRUE;
+ hc->psp.pol_init = TRUE;
+ hc->psp.tor_init = TRUE;
hc->spectral_solution_computed = TRUE;
}
/*
Modified: mc/1D/hc/trunk/hcplates/Makefile
===================================================================
--- mc/1D/hc/trunk/hcplates/Makefile 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hcplates/Makefile 2007-07-12 19:08:51 UTC (rev 7650)
@@ -1,25 +1,66 @@
-CC=gcc
-INCS=hcplates.h
-INCS1=hc_findplate.h
+
+# object file directory
+ODIR = ../../hc/objects/hcplates/$(ARCH)/
+#
+#
+# binary directory
+BDIR = ../../hc/bin/$(ARCH)/
+
+INCS1=hcplates.h
+INCS2=hc_findplate.h
INCS2=hc_ptrot.h
+
+INCS = $(INCS1) $(INCS2) $(INCS3)
+
LIBS=-lm
-LIST = hcplates.c hcplates_init.c zero_arrays.c plates_propagator.c get_loads.c read_unitrots_get_coeff.c integrate_torques.c crlb_vector_plm.c solve_rot.c svdcmp.c svbksb.c crlb_residual.c hcplates_write_output.c
-LIST2 = findplate.c
-LIST3 = ptrot.c hc_ptrot_mem.c hc_ptrot_init.c
-all: hcplates hc_findplate hc_ptrot
+LIST = $(ODIR)/hcplates.o $(ODIR)/hcplates_init.o $(ODIR)/zero_arrays.o \
+ $(ODIR)/plates_propagator.o $(ODIR)/get_loads.o \
+ $(ODIR)/read_unitrots_get_coeff.o $(ODIR)/integrate_torques.o \
+ $(ODIR)/crlb_vector_plm.o $(ODIR)/solve_rot.o $(ODIR)/svdcmp.o \
+ $(ODIR)/svbksb.o $(ODIR)/crlb_residual.o $(ODIR)/hcplates_write_output.o
-hcplates: $(INCS) $(LIST)
- ${CC} $(INCS) -o hcplates $(LIST) $(CFLAGS) $(LIBS)
+LIST2 = $(ODIR)/findplate.o
-hc_findplate: $(INCS1) $(LIST2)
+LIST3 = $(ODIR)/ptrot.o $(ODIR)/hc_ptrot_mem.o $(ODIR)/hc_ptrot_init.o
- gcc -g $(INCS1) -o hc_findplate $(LIST2)
+all: dirs bins
-hc_ptrot: $(INCS2) $(LIST3)
+dirs:
+ if [ ! -s ../../hc/ ]; then\
+ mkdir ../../hc/;\
+ fi;\
+ if [ ! -s ../../hc/objects/ ]; then\
+ mkdir ../../hc/objects;\
+ fi;\
+ if [ ! -s ../../hc/objects/hcplates/ ]; then\
+ mkdir ../../hc/objects/hcplates/;\
+ fi;\
+ if [ ! -s ../../hc/objects/hcplates/$(ARCH)/ ]; then\
+ mkdir ../../hc/objects/hcplates/$(ARCH);\
+ fi;\
+ if [ ! -s ../../hc/bin/ ];then\
+ mkdir ../../hc/bin;\
+ fi;\
+ if [ ! -s ../../hc/bin/$(ARCH) ];then\
+ mkdir ../../hc/bin/$(ARCH);\
+ fi
- gcc -g $(INCS2) -o hc_ptrot $(LIST3)
+bins: $(BDIR)/hcplates $(BDIR)/hc_findplate $(BDIR)/hc_ptrot
+$(BDIR)/hcplates: $(LIST)
+ $(CC) -o $(BDIR)/hcplates $(LIST) $(LIBS)
+
+$(BDIR)/hc_findplate: $(LIST2)
+ $(CC) -o $(BDIR)/hc_findplate $(LIST2)
+
+$(BDIR)/hc_ptrot: $(INCS2) $(LIST3)
+ $(CC) -o $(BDIR)/hc_ptrot $(LIST3)
+
+$(ODIR)/%.o: %.c $(INCS)
+ $(CC) $(CFLAGS) $(INC_FLAGS) $(DEFINES) -c $< -o $(ODIR)/$*.o
+
+
Modified: mc/1D/hc/trunk/hcplates/findplate.c
===================================================================
--- mc/1D/hc/trunk/hcplates/findplate.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hcplates/findplate.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -221,10 +221,8 @@
fprintf(stderr,"options:\n\n");
fprintf(stderr,"-B plate boundary file (%s)\n",
p->boundary_file);
- fprintf(stderr,"-D data file (enes) containg no. plates (first line), then no. points on each plate (rest). Both files needed. \n",
- p->data_filename);
- fprintf(stderr,"-O output file (default: plate_ids.ixz) containg plates id, lat & long. \n",
- p->data_filename);
+ fprintf(stderr,"-D data file (enes) containg no. plates (first line), then no. points on each plate (rest). Both files needed. \n");
+ fprintf(stderr,"-O output file (default: plate_ids.ixz) containg plates id, lat & long\n");
fprintf(stderr,"\n\n");
exit(-1);
}
Modified: mc/1D/hc/trunk/hcplates/hc_ptrot_init.c
===================================================================
--- mc/1D/hc/trunk/hcplates/hc_ptrot_init.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hcplates/hc_ptrot_init.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -96,14 +96,10 @@
fprintf(stderr,"Options:\n\n");
fprintf(stderr,"-B plate boundary file (%s)\n",
p->boundary_file);
- fprintf(stderr,"-D data file (enes) containg no. plates (first line), then no. points on each plate (rest). Both files needed. \n",
- p->data_filename);
- fprintf(stderr,"-P plate id file (default: plate_ids.ixz) containg plates id, lat & long. \n",
- p->data_filename);
- fprintf(stderr,"-O output file (default: unitrots.coeff) containg coefficients for unit rotations on each plate \n",
- p->data_filename);
- fprintf(stderr,"-d degree of spherical harmonics (defaults to 20) \n",
- p->data_filename);
+ fprintf(stderr,"-D data file (enes) containg no. plates (first line), then no. points on each plate (rest). Both files needed.\n");
+ fprintf(stderr,"-P plate id file (default: plate_ids.ixz) containg plates id, lat & long. \n");
+ fprintf(stderr,"-O output file (default: unitrots.coeff) containg coefficients for unit rotations on each plate \n");
+ fprintf(stderr,"-d degree of spherical harmonics (defaults to 20) \n");
fprintf(stderr,"\n\n");
exit(-1);
}
Modified: mc/1D/hc/trunk/hcplates/hcplates_init.c
===================================================================
--- mc/1D/hc/trunk/hcplates/hcplates_init.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hcplates/hcplates_init.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -364,23 +364,23 @@
{ printf("Error loading parameter file. Will use defaults. Ciao. \n"); }
else {
- fscanf(filePtr_P,"%[^\n]%*1[\n]",&toss1); //Skip first 4 lines
+ fscanf(filePtr_P,"%[^\n]%*1[\n]",toss1); //Skip first 4 lines
//fprintf(stderr,"%s \n",toss1);
- fscanf(filePtr_P,"%[^\n]%*1[\n]",&toss1);
+ fscanf(filePtr_P,"%[^\n]%*1[\n]",toss1);
//fprintf(stderr,"%s \n",toss1);
- fscanf(filePtr_P,"%[^\n]%*1[\n]",&toss1);
+ fscanf(filePtr_P,"%[^\n]%*1[\n]",toss1);
//fprintf(stderr,"%s \n",toss1);
- fscanf(filePtr_P,"%[^\n]%*1[\n]",&toss1);
+ fscanf(filePtr_P,"%[^\n]%*1[\n]",toss1);
//fprintf(stderr,"%s \n",toss1);
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
// read no. distinct layers
- fscanf(filePtr_P,"%s %i",&toss1,&nlayer);
+ fscanf(filePtr_P,"%s %i",toss1,&nlayer);
fprintf(stderr," nlayer = %i \n",nlayer);
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
//Scan in viscosity for different layers
- fscanf(filePtr_P,"%s",&toss1);
+ fscanf(filePtr_P,"%s",toss1);
for (i=0; i<=nlayer; i++) {
fscanf(filePtr_P,"%lf",&r2[i]);
}
@@ -391,47 +391,47 @@
}
//fprintf(stderr," %s %lf, %lf, %lf, %lf %lf, %lf\n",toss1,p->r[0],p->r[1],p->r[2],p->r[3],p->r[4],p->r[5]);
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
// Read viscosity
- fscanf(filePtr_P,"%s %lf",&toss1,&p->visc0);
+ fscanf(filePtr_P,"%s %lf",toss1,&p->visc0);
fprintf(stderr," visc0 = %lf \n",p->visc0);
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
// Viscosity layers
- fscanf(filePtr_P,"%s",&toss1);
+ fscanf(filePtr_P,"%s",toss1);
for (i=0; i<nlayer; i++) {
fscanf(filePtr_P,"%lf",&visc2[i]);
fprintf(stderr," %i %lf\n",i,visc2[i]);
}
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
// core radius
- fscanf(filePtr_P,"%s %lf",&toss1,&p->RCore);
+ fscanf(filePtr_P,"%s %lf",toss1,&p->RCore);
fprintf(stderr," Rcore = %lf km \n",p->RCore);
p->RCore *= 1e5; //units not kms
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
// earth radius
- fscanf(filePtr_P,"%s %lf",&toss1,&p->erad);
+ fscanf(filePtr_P,"%s %lf",toss1,&p->erad);
fprintf(stderr," Rearth = %lf km \n",p->erad);
p->erad *= 1e5; //units not kms
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
//No .plates
- fscanf(filePtr_P,"%s %i",&toss1,&p->NPLT);
+ fscanf(filePtr_P,"%s %i",toss1,&p->NPLT);
fprintf(stderr," NPLT = %i \n",p->NPLT);
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
// Plate areas
// Viscosity layers
- fscanf(filePtr_P,"%s",&toss1);
+ fscanf(filePtr_P,"%s",toss1);
for (i=0; i<=p->NPLT; i++) {
fscanf(filePtr_P,"%lf",&p_area[i]);
fprintf(stderr," %i %lf\n",i,p_area[i]);
}
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
//LDIM
- fscanf(filePtr_P,"%s %i",&toss1,&p->LDIM);
+ fscanf(filePtr_P,"%s %i",toss1,&p->LDIM);
fprintf(stderr," LDIM = %i \n",p->LDIM);
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
// Dependents
p->NPDIM = 3 * p->NPLT;
@@ -439,30 +439,30 @@
p->KDIM = (p->LDIM + 2)*(p->LDIM+1)/2;
// Geopgraphic divisor
- fscanf(filePtr_P,"%s %i",&toss1,&p->dd);
+ fscanf(filePtr_P,"%s %i",toss1,&p->dd);
fprintf(stderr," dd = %i \n",p->dd);
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
// time conversion factor
- fscanf(filePtr_P,"%s %lf",&toss1,&p->stoy);
+ fscanf(filePtr_P,"%s %lf",toss1,&p->stoy);
fprintf(stderr," stoy = %lf \n",p->stoy);
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
//Lmax
- fscanf(filePtr_P,"%s %i",&toss1,&p->Lmax);
+ fscanf(filePtr_P,"%s %i",toss1,&p->Lmax);
fprintf(stderr," Lmax = %i \n",p->Lmax);
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
//n
- fscanf(filePtr_P,"%s %i",&toss1,&p->n);
+ fscanf(filePtr_P,"%s %i",toss1,&p->n);
fprintf(stderr," n = %i \n",p->n);
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
//no-slip
- fscanf(filePtr_P,"%s %i",&toss1,&p->iba);
+ fscanf(filePtr_P,"%s %i",toss1,&p->iba);
fprintf(stderr," iba = %i \n",p->iba);
// scaling ratio for velocities at end
//no-slip
- fscanf(filePtr_P," %[^\n]%*1[\n]",&toss1); //This will read blank line and next text line as one
- fscanf(filePtr_P,"%s %lf",&toss1,&p->ratio);
+ fscanf(filePtr_P," %[^\n]%*1[\n]",toss1); //This will read blank line and next text line as one
+ fscanf(filePtr_P,"%s %lf",toss1,&p->ratio);
fprintf(stderr," ratio = %lf \n",p->ratio);
//DONE
Modified: mc/1D/hc/trunk/hcplates/ptrot.c
===================================================================
--- mc/1D/hc/trunk/hcplates/ptrot.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hcplates/ptrot.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -253,7 +253,8 @@
double dvphixp,dvphixm,dvphiyp,dvphiym,dvphizp,dvphizm;
double term1x,term1y,term1z,term2xp,term2yp,term2zp,term2xm,term2ym,term2zm;
double term3xp,term3yp,term3zp,term3xm,term3ym,term3zm,theta,avephi,x,facx;
-
+ /* added this TWB */
+ toss4 = 360;
/*Opening up boundary file */
strncpy(file4,p->plateids_file,HC_CHAR_LENGTH);
filePtr4 = fopen(file4, "r");
Modified: mc/1D/hc/trunk/hcplates/run_example.mk
===================================================================
--- mc/1D/hc/trunk/hcplates/run_example.mk 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/hcplates/run_example.mk 2007-07-12 19:08:51 UTC (rev 7650)
@@ -6,4 +6,6 @@
hcplates -L point.j -T plates_ids.ixz -U unitrots_new.coeffs -Op mypoles.out -Ov myvels.out -P parameter_file.default
+pmyvels_simple
+
Modified: mc/1D/hc/trunk/main.c
===================================================================
--- mc/1D/hc/trunk/main.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/main.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -127,14 +127,14 @@
*/
/*
- make room for the spectral solution
+ make room for the spectral solution on irregular grid
*/
sh_allocate_and_init(&sol_spectral,nsol,lmax,model->sh_type,1,
- p->verbose);
+ p->verbose,FALSE);
if(p->compute_geoid)
/* make room for geoid solution */
sh_allocate_and_init(&geoid,1,model->dens_anom[0].lmax,
- model->sh_type,0,p->verbose);
+ model->sh_type,0,p->verbose,FALSE);
/*
solve poloidal and toroidal part and sum
@@ -214,7 +214,7 @@
free(sol_spatial);
if(p->verbose)
fprintf(stderr,"%s: done\n",argv[0]);
-
+ hc_struc_free(&model);
return 0;
}
Modified: mc/1D/hc/trunk/pow.sh_ana
===================================================================
--- mc/1D/hc/trunk/pow.sh_ana 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/pow.sh_ana 2007-07-12 19:08:51 UTC (rev 7650)
@@ -1,32 +1,32 @@
- 0 1.0000005e+00
- 1 1.6666666e+00
- 2 1.8000000e+00
- 3 1.8571428e+00
- 4 4.5955911e-14
- 5 1.2534232e-13
- 6 7.0812700e-14
- 7 7.8964782e-14
- 8 5.4092612e-14
- 9 1.3155515e-13
- 10 8.5687845e-14
- 11 8.9160592e-14
- 12 3.8386463e-14
- 13 8.9618200e-14
- 14 6.0837667e-14
- 15 8.9549502e-14
- 16 3.2622397e-14
- 17 6.5555946e-14
- 18 3.6024145e-14
- 19 6.9451592e-14
- 20 2.7890501e-14
- 21 5.0478410e-14
- 22 3.5044345e-14
- 23 8.8381484e-14
- 24 4.3674510e-14
- 25 3.7003126e-14
- 26 3.4498568e-14
- 27 7.7315819e-14
- 28 5.0666421e-14
- 29 4.6164197e-14
- 30 4.2085591e-14
- 31 6.4966384e-14
+ 0 0.0000000e+00 0.0000000e+00
+ 1 1.0039276e+00 1.0039278e+00
+ 2 1.0103023e+00 1.0103025e+00
+ 3 9.9611825e-01 9.9611843e-01
+ 4 1.1189895e-05 1.1189598e-05
+ 5 1.3912736e-06 1.3914931e-06
+ 6 3.6852563e-07 3.6840964e-07
+ 7 1.4058963e-07 1.4065725e-07
+ 8 6.4872495e-08 6.4843256e-08
+ 9 3.4803097e-08 3.4802554e-08
+ 10 2.0175317e-08 2.0168704e-08
+ 11 1.2711423e-08 1.2716271e-08
+ 12 8.3854665e-09 8.3869294e-09
+ 13 5.8351342e-09 5.8357870e-09
+ 14 4.1510799e-09 4.1510146e-09
+ 15 3.1080256e-09 3.1102638e-09
+ 16 2.3411473e-09 2.3384403e-09
+ 17 1.8074796e-09 1.8086181e-09
+ 18 1.4325702e-09 1.4303744e-09
+ 19 1.1559325e-09 1.1544544e-09
+ 20 9.3577301e-10 9.3553199e-10
+ 21 7.8341189e-10 7.8380458e-10
+ 22 6.5448347e-10 6.5526079e-10
+ 23 5.5339594e-10 5.5376032e-10
+ 24 4.7935256e-10 4.7784049e-10
+ 25 4.1667336e-10 4.1782511e-10
+ 26 3.6498404e-10 3.6626421e-10
+ 27 3.2285002e-10 3.2188649e-10
+ 28 2.9063027e-10 2.9093525e-10
+ 29 2.6214117e-10 2.6328723e-10
+ 30 2.3942193e-10 2.3912286e-10
+ 31 2.2615471e-10 2.2604703e-10
Modified: mc/1D/hc/trunk/pow.shana
===================================================================
--- mc/1D/hc/trunk/pow.shana 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/pow.shana 2007-07-12 19:08:51 UTC (rev 7650)
@@ -1,32 +1,32 @@
-0 9.99958600428490e-01
-1 9.99967660778602e-01
-2 9.99981440424864e-01
-3 9.99979000772176e-01
-4 4.20608703170892e-10
-5 7.71255339458584e-10
-6 4.21497056283272e-10
-7 7.71715943934347e-10
-8 4.21676338450171e-10
-9 7.71913040724380e-10
-10 4.21784309298177e-10
-11 7.72075570761280e-10
-12 4.21875455706380e-10
-13 7.72238808731481e-10
-14 4.21968554816294e-10
-15 7.72416459219568e-10
-16 4.22071288085273e-10
-17 7.72613821070832e-10
-18 4.22188997253500e-10
-19 7.72835594627370e-10
-20 4.22315060422342e-10
-21 7.73080383162030e-10
-22 4.22456780035087e-10
-23 7.73348357567720e-10
-24 4.22609744865580e-10
-25 7.73636353330838e-10
-26 4.22771325874776e-10
-27 7.73948786617862e-10
-28 4.22951878724432e-10
-29 7.74289603722983e-10
-30 4.23142562338938e-10
-31 7.74648627700129e-10
+0 0.00000000000000e+00 0.00000000000000e+00
+1 1.00001549851264e+00 1.00001554178822e+00
+2 1.00001061977234e+00 1.00001061865053e+00
+3 1.00001292930198e+00 1.00001294221066e+00
+4 2.98584189355089e-11 3.00499293609165e-11
+5 2.96210939023542e-11 2.96959378792701e-11
+6 2.97470149011265e-11 2.98107019757459e-11
+7 2.96968479310751e-11 2.97366325374235e-11
+8 2.97405582599865e-11 2.97740142538467e-11
+9 2.97229644755101e-11 2.97436763405639e-11
+10 2.97427336140716e-11 2.97630885605110e-11
+11 2.97342663571140e-11 2.97505370056352e-11
+12 2.97468125222950e-11 2.97605917946227e-11
+13 2.97434659496090e-11 2.97528816496266e-11
+14 2.97523460897165e-11 2.97627319621183e-11
+15 2.97511102392692e-11 2.97597571690397e-11
+16 2.97585384159843e-11 2.97661210544403e-11
+17 2.97595708459170e-11 2.97658854760279e-11
+18 2.97659486431519e-11 2.97716585105063e-11
+19 2.97685377364341e-11 2.97722967135828e-11
+20 2.97739066236587e-11 2.97782682503121e-11
+21 2.97759259172495e-11 2.97805563713286e-11
+22 2.97824264986484e-11 2.97865495622269e-11
+23 2.97860435904979e-11 2.97891843206281e-11
+24 2.97917917193661e-11 2.97952504422914e-11
+25 2.97957862387590e-11 2.97992031252785e-11
+26 2.98016405682321e-11 2.98039173370977e-11
+27 2.98063481713644e-11 2.98085704276222e-11
+28 2.98127253998186e-11 2.98153085308753e-11
+29 2.98180935434980e-11 2.98200348075614e-11
+30 2.98250181704652e-11 2.98265130279608e-11
+31 2.98300931845548e-11 2.98323483303453e-11
Modified: mc/1D/hc/trunk/rick_sh_c.c
===================================================================
--- mc/1D/hc/trunk/rick_sh_c.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/rick_sh_c.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -6,8 +6,8 @@
// ivec is set to 1
// */
-void rick_compute_allplm(int lmax,int ivec,double *plm,
- double *dplm, struct rick_module *rick)
+void rick_compute_allplm(int lmax,int ivec,float *plm,
+ float *dplm, struct rick_module *rick)
{
int i,os;
@@ -17,19 +17,36 @@
exit(-1);
}
os=0; /* changed this to 0 TWB */
- for (i=0;i<rick->nlat;i++) { /*changed from 1->nlat to 0->nlat-1 - need change in plmbar1 also */
+ for (i=0;i < rick->nlat;i++) { /*changed from 1->nlat to 0->nlat-1 - need change in plmbar1 also */
rick_plmbar1((plm+os),(dplm+os),ivec,lmax,rick->gauss_z[i],rick); /*note change in gauss_z[i] */
os += rick->lmsize;
}
}
+/* // compute Legendre function (l,m) evaluated on npoints points in
+// theta array their derviatives with respect to theta, if ivec is set
+// to 1 */
+void rick_compute_allplm_irreg(int lmax,int ivec,SH_RICK_PREC *plm,
+ SH_RICK_PREC *dplm, struct rick_module *rick,
+ float *theta, int ntheta)
+{
+ int i,os;
+
+ os=0; /* changed this to 0 TWB */
+ for (i=0;i < ntheta;i++) { /*changed from 1->nlat to 0->nlat-1 - need change in plmbar1 also */
+ rick_plmbar1((plm+os),(dplm+os),ivec,lmax,cos(theta[i]),rick); /*note change in gauss_z[i] */
+ os += rick->lmsize;
+ }
+}
+
+
/* //
// detemine the colatidude and longitude of PIxel index
// where index goes from 0 ... nlon * nlat-1
// */
-void rick_pix2ang(int index, int lmax, double *theta,
- double *phi, struct rick_module *rick) {
+void rick_pix2ang(int index, int lmax, SH_RICK_PREC *theta,
+ SH_RICK_PREC *phi, struct rick_module *rick) {
int i,j;
if(!rick->initialized){
fprintf(stderr,"rick_pix2ang: error: rick module not initialized\n");
@@ -79,7 +96,7 @@
//
// input: lmax,ivec
// local
- double *plm,*dplm;
+ SH_RICK_PREC *plm,*dplm;
if(!rick->initialized){
fprintf(stderr,"rick_shc2d: error: initialize first\n");
exit(-1);
@@ -91,9 +108,9 @@
exit(-1);
}
/* allocate memory */
- hc_dvecalloc(&plm,rick->nlat*rick->lmsize,"rick_shc2d: mem 1");
+ hc_svecalloc(&plm,rick->nlat*rick->lmsize,"rick_shc2d: mem 1");
if(ivec)
- hc_dvecalloc(&dplm,rick->nlat*rick->lmsize,"rick_shc2d: mem 2");
+ hc_svecalloc(&dplm,rick->nlat*rick->lmsize,"rick_shc2d: mem 2");
//
// compute the Plm first
rick_compute_allplm(lmax,ivec,plm,dplm,rick);
@@ -107,14 +124,47 @@
if(ivec)
free(dplm);
}
+/*
+converts on irregular basis with locations cos(theta)[], phi[] long
+ */
+void rick_shc2d_irreg(SH_RICK_PREC *cslm,SH_RICK_PREC *dslm,
+ int lmax,int ivec,
+ SH_RICK_PREC *rdatax,SH_RICK_PREC *rdatay,
+ struct rick_module *rick,float *theta, int ntheta,
+ float *phi,int nphi, hc_boolean save_sincos_fac)
+{
+ SH_RICK_PREC *plm,*dplm;
+ if(!rick->initialized){
+ fprintf(stderr,"rick_shc2d: error: initialize first\n");
+ exit(-1);
+ }
+ /* allocate memory */
+ hc_svecalloc(&plm,ntheta*rick->lmsize,"rick_shc2d_irreg: mem 1");
+ if(ivec)
+ hc_svecalloc(&dplm,ntheta*rick->lmsize,"rick_shc2d_irreg: mem 2");
+ //
+ // compute the Plm first
+ rick_compute_allplm_irreg(lmax,ivec,plm,dplm,rick,theta,ntheta);
+ //
+ // call the precomputed subroutine
+ rick_shc2d_pre_irreg(cslm,dslm,lmax,plm,dplm,ivec,rdatax,rdatay,rick,theta,ntheta,
+ phi,nphi,save_sincos_fac);
+ /* free legendre functions if not needed */
+ free(plm);
+ if(ivec)
+ free(dplm);
+}
+
+
+
/* //
// the actual routine to go from spectral to spatial: added structure rick
// */
void rick_shc2d_pre(SH_RICK_PREC *cslm,SH_RICK_PREC *dslm,
- int lmax,double *plm, double *dplm,
+ int lmax,SH_RICK_PREC *plm, SH_RICK_PREC *dplm,
int ivec,float *rdatax,float *rdatay,
struct rick_module *rick)
{
@@ -122,8 +172,8 @@
// Legendre functions are precomputed
// */
SH_RICK_PREC *valuex, *valuey;
- double dpdt,dpdp;
- static int negunity = -1;
+ SH_RICK_PREC dpdt,dpdp;
+ static int negunity = -1; /* an actual constant */
int i,j,m,m2,j2,ios1,l,lmaxp1,lmaxp1t2,oplm,nlon2,lm1;
if(!rick->initialized){
fprintf(stderr,"rick_shc2d_pre: error: initialize modules first\n");
@@ -172,9 +222,9 @@
//fprintf(stderr,"%11i %11i %11g %11g\n",l,m,cslm[j2],cslm[j2+1]);
/* add up contributions from all l,m */
valuex[m2] += /* cos term */
- plm[oplm+j] * (double)cslm[j2]; /* A coeff */
+ plm[oplm+j] * (SH_RICK_PREC)cslm[j2]; /* A coeff */
valuex[m2+1] += /* sin term */
- plm[oplm+j] * (double)cslm[j2+1]; /* B coeff */
+ plm[oplm+j] * (SH_RICK_PREC)cslm[j2+1]; /* B coeff */
}
/* compute inverse FFT */
@@ -210,9 +260,9 @@
m2 = 2*m;
/* derivative factors */
lm1 = l - 1;
- dpdt = dplm[oplm+j] * (double)rick->ell_factor[lm1]; /* d_theta(P_lm) factor */
- dpdp = ((double)m) * plm[oplm+j]/ (double)rick->sin_theta[i];
- dpdp *= (double)rick->ell_factor[lm1]; /* d_phi (P_lm) factor */
+ dpdt = dplm[oplm+j] * (SH_RICK_PREC)rick->ell_factor[lm1]; /* d_theta(P_lm) factor */
+ dpdp = ((SH_RICK_PREC)m) * plm[oplm+j]/ (SH_RICK_PREC)rick->sin_theta[i];
+ dpdp *= (SH_RICK_PREC)rick->ell_factor[lm1]; /* d_phi (P_lm) factor */
/* add up contributions from all l,m
u_theta
@@ -257,6 +307,140 @@
}
+/* //
+
+irregular version
+
+// */
+void rick_shc2d_pre_irreg(SH_RICK_PREC *cslm,SH_RICK_PREC *dslm,
+ int lmax,SH_RICK_PREC *plm, SH_RICK_PREC *dplm,
+ int ivec,float *rdatax,float *rdatay,
+ struct rick_module *rick, float *theta, int ntheta,
+ float *phi,int nphi,my_boolean save_sincos_fac)
+{
+ /* //
+ // Legendre functions are precomputed
+ // */
+ double dpdt,dpdp,mphi,sin_theta;
+ float *loc_plma=NULL,*loc_plmb=NULL;
+ int i,j,k,k2,m,ios1,ios2,ios3,l,lmaxp1,lm1,idata;
+ if(!rick->initialized){
+ fprintf(stderr,"rick_shc2d_pre_irreg: error: initialize modules first\n");
+ exit(-1);
+ }
+
+ lmaxp1 = lmax + 1; // this is nlat
+
+ if(ivec){
+ if(!rick->vector_sh_fac_init){
+ fprintf(stderr,"rick_shc2d_pre_irreg: error: vector harmonics factors not initialized\n");
+ exit(-1);
+ }
+ }
+ if((lmax+1)*(lmax+2)/2 > rick->lmsize){
+ fprintf(stderr,"rick_shc2d_pre_irreg: error: lmax %i out of bounds\n",lmax);
+ exit(-1);
+ }
+
+ /*
+ compute sin/cos factors
+ */
+ hc_svecrealloc(&rick->sfac,nphi*lmaxp1,"rick_shc2d_pre_irreg");
+ hc_svecrealloc(&rick->cfac,nphi*lmaxp1,"rick_shc2d_pre_irreg");
+ for(ios1=i=0;i < nphi;i++){
+ for(m=0;m <= lmax;m++,ios1++){
+ mphi = (double)m*(double)phi[i];
+ rick->sfac[ios1] = (float)sin(mphi);
+ rick->cfac[ios1] = (float)cos(mphi);
+ }
+ }
+ if(ivec == 0){
+ /*
+
+ scalar
+
+ */
+ hc_svecrealloc(&loc_plma,ntheta*rick->lmsize,"rick_shc2d_pre_irreg 3");
+ hc_svecrealloc(&loc_plmb,ntheta*rick->lmsize,"rick_shc2d_pre_irreg 4");
+ for(i=ios1=0;i < ntheta;i++){
+ for(k=k2=0;k < rick->lmsize;k++,k2+=2,ios1++){
+ loc_plma[ios1] = cslm[k2 ] * plm[ios1];
+ loc_plmb[ios1] = cslm[k2+1] * plm[ios1];
+ }
+ }
+ for (idata=i=ios2=0;i < ntheta; i++,ios2 += rick->lmsize) { /* theta
+ loop */
+ for(ios3=j=0;j < nphi;j++,idata++,ios3 += lmaxp1){ /* phi loop */
+
+ /* m = 0 , l = 0*/
+ l=0;m=0;
+ rdatax[idata] = loc_plma[ios2];
+ for (k=1; k < rick->lmsize; k++) {
+ m++;
+ if (m > l) {
+ m=0;l++;
+ }
+ rdatax[idata] += loc_plma[ios2+k] * rick->cfac[ios3+m];
+ rdatax[idata] += loc_plmb[ios2+k] * rick->sfac[ios3+m];
+ }
+ }
+ }
+
+ free(loc_plma);free(loc_plmb);
+ /* end scalar part */
+ } else {
+ /*
+ vector harmonics
+ */
+ for (idata=i=ios2=0;i < ntheta; i++,ios2 += rick->lmsize) { /* theta
+ loop */
+ sin_theta = sin(theta[i]);
+ for(ios3=j=0;j < nphi;j++,idata++,ios3 += lmaxp1){ /* phi loop */
+
+ rdatax[idata] = rdatay[idata] = 0.0;
+
+ l = 0;m = -1;
+ /* start at l = 1 */
+ for (k=1,k2=2; k < rick->lmsize; k++,k2+=2) {
+ /* loop through l,m */
+ m++;
+ if (m > l) {
+ m=0;l++;
+ }
+ lm1 = l - 1;
+ dpdt = dplm[ios2+k] * rick->ell_factor[lm1]; /* d_theta(P_lm) factor */
+ dpdp = ((SH_RICK_PREC)m) * plm[ios2+k]/ sin_theta;
+ dpdp *= (SH_RICK_PREC)rick->ell_factor[lm1]; /* d_phi (P_lm) factor */
+
+ /*
+
+ add up contributions from all l,m
+
+ u_theta
+
+ */
+ rdatax[i] += /* cos term */
+ (cslm[k2] * dpdt + dslm[k2+1] * dpdp) * rick->cfac[ios3+m];
+ rdatax[i] += /* sin term */
+ (cslm[k2+1] * dpdt - dslm[k2] * dpdp) * rick->sfac[ios3+m];
+ /*
+ u_phi
+ */
+ rdatay[i] += // cos term
+ (cslm[k2+1] * dpdp - dslm[k2] * dpdt) * rick->cfac[ios3+m];
+ rdatay[i] += // sin term
+ (- cslm[k2] * dpdp - dslm[k2+1] * dpdt) * rick->sfac[ios3+m];
+ }
+ }
+ }
+ }
+
+ if(!save_sincos_fac){
+ hc_svecrealloc(&rick->sfac,1,"");
+ hc_svecrealloc(&rick->cfac,1,"");
+ }
+}
+
void rick_shd2c(SH_RICK_PREC *rdatax,SH_RICK_PREC *rdatay,
int lmax,int ivec,SH_RICK_PREC *cslm,
SH_RICK_PREC *dslm,struct rick_module *rick)
@@ -307,10 +491,10 @@
//
//
// local
- double *plm,*dplm;
+ SH_RICK_PREC *plm,*dplm;
/* allocate memory */
- hc_dvecalloc(&plm,rick->nlat*rick->lmsize,"rick_shd2c: mem 1");
- hc_dvecalloc(&dplm,rick->nlat*rick->lmsize,"rick_shd2c: mem 2");
+ hc_svecalloc(&plm,rick->nlat*rick->lmsize,"rick_shd2c: mem 1");
+ hc_svecalloc(&dplm,rick->nlat*rick->lmsize,"rick_shd2c: mem 2");
// check
if(!rick->initialized){
fprintf(stderr,"rick_shd2c: error: initialize first\n");
@@ -336,14 +520,14 @@
// for comments, see above
//
void rick_shd2c_pre(SH_RICK_PREC *rdatax,SH_RICK_PREC *rdatay,
- int lmax,double *plm,double *dplm,int ivec,
+ int lmax,SH_RICK_PREC *plm,SH_RICK_PREC *dplm,int ivec,
SH_RICK_PREC *cslm,SH_RICK_PREC *dslm,
struct rick_module *rick)
{
// local
SH_RICK_PREC *valuex, *valuey;
- double dfact,dpdt,dpdp;
- static int unity = 1;
+ SH_RICK_PREC dfact,dpdt,dpdp;
+ static int unity = 1; /* constant */
//
int lmaxp1,lmaxp1t2,i,j,l,m,ios1,m2,j2,oplm,nlon2,lm1;
// check
@@ -414,9 +598,9 @@
}
// we incorporate the Gauss integration weight and Plm factors here
if (m == 0) {
- dfact = ((double)rick->gauss_w[i] * plm[oplm+j])/2.0;
+ dfact = ((SH_RICK_PREC)rick->gauss_w[i] * plm[oplm+j])/2.0;
}else{
- dfact = ((double)rick->gauss_w[i] * plm[oplm+j])/4.0;
+ dfact = ((SH_RICK_PREC)rick->gauss_w[i] * plm[oplm+j])/4.0;
}
m2 = m * 2;
cslm[j2] += valuex[m2] * dfact; // A coefficient
@@ -463,7 +647,7 @@
// d_theta(P_lm) factor
dpdt = dplm[oplm+j];
// d_phi (P_lm) factor
- dpdp = ((double)m) * plm[oplm+j]/(double)rick->sin_theta[i];
+ dpdp = ((SH_RICK_PREC)m) * plm[oplm+j]/(SH_RICK_PREC)rick->sin_theta[i];
//
m2 = m * 2;
// print *,m,l,dpdt*dfact,dpdp*dfact
@@ -492,8 +676,11 @@
//
// if ivec == 1, will initialize for velocities/polarizations
//
+/* if regular is set, will not initialize the Gauss quadrature
+ points */
void rick_init(int lmax,int ivec,int *npoints,int *nplm,
- int *tnplm, struct rick_module *rick)
+ int *tnplm, struct rick_module *rick,
+ hc_boolean regular)
{
// input: lmax,ivec
@@ -501,8 +688,8 @@
// local
SH_RICK_PREC xtemp;
int i,l;
- static int old_lmax,old_ivec,old_npoints,old_nplm,old_tnplm;
+
if(!rick->was_called){
if(lmax == 0){
fprintf(stderr,"rick_init: error: lmax is zero: %i\n",lmax);
@@ -528,7 +715,7 @@
// number of points in one layer
//
*npoints = rick->nlat * rick->nlon;
- old_npoints = *npoints;
+ rick->old_npoints = *npoints;
//
//
// for coordinate computations
@@ -544,47 +731,52 @@
// size of the Plm array
*nplm = rick->lmsize * rick->nlat;
*tnplm = *nplm * (1+ivec); // for all layers
- old_tnplm = *tnplm;
- old_nplm = *nplm;
+ rick->old_tnplm = *tnplm;
+ rick->old_nplm = *nplm;
//
// initialize the Gauss points, at which the latitudes are
// evaluated
//
- rick_vecalloc(&rick->gauss_z,rick->nlat,"rick_init 1");
- rick_vecalloc(&rick->gauss_w,rick->nlat,"rick_init 2");
- rick_vecalloc(&rick->gauss_theta,rick->nlat,"rick_init 3");
- /*
- gauss weighting
- */
- rick_gauleg(-1.0,1.0,rick->gauss_z,rick->gauss_w,rick->nlat);
+ if(!regular){
+ rick_vecalloc(&rick->gauss_z,rick->nlat,"rick_init 1");
+ rick_vecalloc(&rick->gauss_w,rick->nlat,"rick_init 2");
+ rick_vecalloc(&rick->gauss_theta,rick->nlat,"rick_init 3");
+ /*
+ gauss weighting
+ */
+ rick_gauleg(-1.0,1.0,rick->gauss_z,rick->gauss_w,rick->nlat);
+ //
+ // theta values of the Gauss quadrature points
+ //
+ for(i=0;i < rick->nlat;i++)
+ rick->gauss_theta[i] = acos(rick->gauss_z[i]);
+ }
//
- // theta values of the Gauss quadrature points
- //
- for(i=0;i < rick->nlat;i++)
- rick->gauss_theta[i] = acos(rick->gauss_z[i]);
- //
// those will be used by plmbar to store some of the factors
//
- hc_dvecalloc(&rick->plm_f1,rick->lmsize,"rick_init 4");
- hc_dvecalloc(&rick->plm_f2,rick->lmsize,"rick_init 5");
- hc_dvecalloc(&rick->plm_fac1,rick->lmsize,"rick_init 6");
- hc_dvecalloc(&rick->plm_fac2,rick->lmsize,"rick_init 7");
- hc_dvecalloc(&rick->plm_srt,rick->nlon,"rick_init 8");
+ hc_svecalloc(&rick->plm_f1,rick->lmsize,"rick_init 4");
+ hc_svecalloc(&rick->plm_f2,rick->lmsize,"rick_init 5");
+ hc_svecalloc(&rick->plm_fac1,rick->lmsize,"rick_init 6");
+ hc_svecalloc(&rick->plm_fac2,rick->lmsize,"rick_init 7");
+ hc_svecalloc(&rick->plm_srt,rick->nlon,"rick_init 8");
if(ivec){
//
// additional arrays for vector spherical harmonics
// (perform the computations in SH_RICK_PREC precision)
//
rick_vecalloc(&rick->ell_factor,rick->nlat,"rick init 9");
- rick_vecalloc(&rick->sin_theta,rick->nlat,"rick init 9");
+ if(!regular)
+ rick_vecalloc(&rick->sin_theta,rick->nlat,"rick init 9");
// 1/(l(l+1)) factors
for(i=0,l=1;i < rick->nlat;i++,l++){
// no l=0 term, obviously
// go from l=1 to l=lmax+1
rick->ell_factor[i] = 1.0/sqrt((SH_RICK_PREC)(l*(l+1)));
- rick->sin_theta[i] = sqrt((1.0 - rick->gauss_z[i])*
- (1.0+rick->gauss_z[i]));
+ if(!regular){
+ rick->sin_theta[i] = sqrt((1.0 - rick->gauss_z[i])*
+ (1.0+rick->gauss_z[i]));
+ }
rick->vector_sh_fac_init = TRUE;
}
}else{
@@ -599,23 +791,26 @@
/*
save initial call lmax and ivec settings
*/
- old_lmax = lmax;
- old_ivec = ivec;
+ rick->old_lmax = lmax;
+ rick->old_ivec = ivec;
+ /* for irregular expansions */
+ rick->cfac = rick->sfac = NULL;
+
/* end initial branch */
}else{
- if(lmax != old_lmax){
+ if(lmax != rick->old_lmax){
fprintf(stderr,"rick_init: error: was init with lmax %i, now: %i. (ivec: %i, now: %i)\n",
- old_lmax,lmax,old_ivec,ivec);
+ rick->old_lmax,lmax,rick->old_ivec,ivec);
exit(-1);
}
- if(ivec > old_ivec){
- fprintf(stderr,"rick_init: error: original ivec %i, now %i\n",old_ivec,ivec);
+ if(ivec > rick->old_ivec){
+ fprintf(stderr,"rick_init: error: original ivec %i, now %i\n",rick->old_ivec,ivec);
exit(-1);
}
- *npoints = old_npoints;
- *nplm = old_nplm;
- *tnplm = old_tnplm;
+ *npoints = rick->old_npoints;
+ *nplm = rick->old_nplm;
+ *tnplm = rick->old_tnplm;
}
} /* end rick init */
@@ -640,8 +835,9 @@
if(ivec){
free(rick->ell_factor);free(rick->sin_theta);
}
+
}
-void rick_plmbar1(double *p,double *dp,int ivec,int lmax,
+void rick_plmbar1(SH_RICK_PREC *p,SH_RICK_PREC *dp,int ivec,int lmax,
SH_RICK_PREC z, struct rick_module *rick)
{
//
@@ -687,7 +883,6 @@
double plm,pm1,pm2,pmm,sintsq,fnum,fden;
//
int i,l,m,k,kstart,l2,mstop,lmaxm1;
- static int old_lmax,old_ivec;
if(!rick->initialized){
fprintf(stderr,"rick_plmbar1: error: module not initialized, call rick_init first\n");
exit(-1);
@@ -709,7 +904,7 @@
// first call, set up some factors. the arrays were allocated in rick_init
//
for(k=0,i=1;k < rick->nlon;k++,i++){
- rick->plm_srt[k] = sqrt((double)(i));
+ rick->plm_srt[k] = sqrt((SH_RICK_PREC)(i));
}
// initialize plm factors
for(i=0;i < rick->lmsize;i++){
@@ -758,8 +953,8 @@
k++;
}
} /* end ivec==1 */
- old_lmax = lmax;
- old_ivec = ivec;
+ rick->old_lmax = lmax;
+ rick->old_ivec = ivec;
rick->computed_legendre = TRUE;
/*
end first call
@@ -769,13 +964,13 @@
subsequent call
*/
// test if lmax has changed
- if(lmax != old_lmax){
- fprintf(stderr,"rick_plmbar1: error: factors were computed for lmax %in",old_lmax);
+ if(lmax != rick->old_lmax){
+ fprintf(stderr,"rick_plmbar1: error: factors were computed for lmax %in",rick->old_lmax);
fprintf(stderr,"rick_plmbar1: error: now, lmax is %i\n",lmax);
exit(-1);
}
- if(ivec > old_ivec){
- fprintf(stderr,"rick_plmbar1: error: init with %i, now ivec %i\n",old_ivec,ivec);
+ if(ivec > rick->old_ivec){
+ fprintf(stderr,"rick_plmbar1: error: init with %i, now ivec %i\n",rick->old_ivec,ivec);
exit(-1);
}
}
Modified: mc/1D/hc/trunk/sh.h
===================================================================
--- mc/1D/hc/trunk/sh.h 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/sh.h 2007-07-12 19:08:51 UTC (rev 7650)
@@ -65,11 +65,18 @@
/*
number of entries for the Legendre function array
*/
- int n_plm,tn_plm;
+ int n_plm,tn_plm,tn_plm_irr;
/*
number of points in each layer the spatial domain
*/
int npoints;
+ /*
+
+ */
+ hc_boolean plm_computed,plm_computed_irr;
+ int old_lmax,old_ivec,old_tnplm,old_tnplm_irr,old_lmax_irr,old_ivec_irr;
+
+
/*
holds the coefficients:
@@ -115,9 +122,9 @@
/* expansions */
struct sh_lms *exp;
int nexp; /* number of expansions */
- double *plm; /* precomputed Legendre
+ SH_RICK_PREC *plm; /* precomputed Legendre
functions */
-
+
/*
spatial data points
*/
Modified: mc/1D/hc/trunk/sh_ana.c
===================================================================
--- mc/1D/hc/trunk/sh_ana.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/sh_ana.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -20,6 +20,18 @@
where data.lonlatz has the data at the required locations as put out by sh_ana -lmax
+
+for scalars, the input format is
+
+lon[deg] lat[deg] scalar
+
+for vectors
+
+lon[deg] lat[deg] v_theta v_phi
+
+where theta and phi are the vector components in South and East direction, respectively.
+
+
$Id: sh_ana.c,v 1.6 2006/01/22 01:11:34 becker Exp $
*/
@@ -31,7 +43,7 @@
hc_boolean verbose = TRUE, use_3d = FALSE, short_format = FALSE,
binary = FALSE, print_spatial_base = FALSE;
float *data, zlabel = 0,*flt_dummy;
- double *dbl_dummy;
+ SH_RICK_PREC *dummy;
HC_PREC fac[3] = {1.,1.,1.};
/*
@@ -75,7 +87,7 @@
*/
shps = (ivec)?(2):(1);
/* intialize expansion first */
- sh_allocate_and_init(&exp,shps*nset,lmax,type,ivec,verbose);
+ sh_allocate_and_init(&exp,shps*nset,lmax,type,ivec,verbose,FALSE);
/* make room for data */
hc_svecalloc(&data,shps * exp->npoints,"sh_ana");
if(print_spatial_base){
@@ -92,7 +104,7 @@
/*
perform spherical harmonic expansion
*/
- sh_compute_spectral(data,ivec,FALSE,&dbl_dummy,
+ sh_compute_spectral(data,ivec,FALSE,&dummy,
exp,verbose);
/* print parameters of expansion */
sh_print_parameters_to_file(exp,shps,ilayer,nset,zlabel,
Modified: mc/1D/hc/trunk/sh_exp.c
===================================================================
--- mc/1D/hc/trunk/sh_exp.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/sh_exp.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -16,7 +16,8 @@
*/
void sh_allocate_and_init(struct sh_lms **exp, int n,int lmax,
- int type, int ivec, hc_boolean verbose)
+ int type, int ivec, hc_boolean verbose,
+ hc_boolean regular)
{
int i;
/* init as zeroes! */
@@ -24,7 +25,7 @@
if(!(*exp))
HC_MEMERROR("sh_allocate_and_init");
for(i=0;i<n;i++){
- sh_init_expansion((*exp+i),lmax,type,ivec,verbose);
+ sh_init_expansion((*exp+i),lmax,type,ivec,verbose,regular);
}
}
@@ -39,10 +40,11 @@
coefficients are initialized as zero
+if regular is set, will not use Gauss points
*/
void sh_init_expansion(struct sh_lms *exp, int lmax, int type,
- int ivec, hc_boolean verbose)
+ int ivec, hc_boolean verbose,hc_boolean regular)
{
/*
initialize logic flags
@@ -75,7 +77,12 @@
exp->lmsmall2 = (exp->lmaxp1)*(exp->lmaxp1+1); /* for A and B */
/*
+
+ */
+ exp->plm_computed = FALSE;
+ /*
+
allocate the spectral (coefficients) storage and initialize possibly
other arrays
@@ -84,6 +91,8 @@
#ifdef HC_USE_HEALPIX
case SH_HEALPIX: /* SH_HEALPIX part */
+ if(regular)
+ HC_ERROR("regular init not implemented for healpix");
/*
get single precision complex array which holds A and B
@@ -123,16 +132,17 @@
*/
#ifdef NO_RICK_FORTRAN
rick_init(exp->lmax,ivec,&exp->npoints,
- &exp->n_plm,&exp->tn_plm,&exp->rick);
+ &exp->n_plm,&exp->tn_plm,&exp->rick,regular);
#else
/* f90 version */
rick_f90_init(&exp->lmax,&ivec,&exp->npoints,
- &exp->n_plm,&exp->tn_plm);
+ &exp->n_plm,&exp->tn_plm,regular);
#endif
break;
#ifdef HC_USE_SPHEREPACK
case SH_SPHEREPACK_GAUSS:
case SH_SPHEREPACK_EVEN:
+ HC_ERROR("init_expansion","Spherepack not implemented");
/*
make room for coefficients
*/
@@ -169,7 +179,9 @@
#ifdef HC_USE_SPHEREPACK
case SH_SPHEREPACK_GAUSS:
case SH_SPHEREPACK_EVEN:
+ HC_ERROR("free_expansion","Spherepack not implemented");
break;
+
#endif
default:
sh_exp_type_error("sh_free_expansion",(exp+i));
@@ -198,7 +210,7 @@
#ifdef HC_USE_SPHEREPACK
case SH_SPHEREPACK_GAUSS:
case SH_SPHEREPACK_EVEN:
-
+ HC_ERROR("clear_alm","Spherepack not implemented");
break;
#endif
default:
@@ -648,7 +660,7 @@
my_boolean use_3d,
int shps, float *data, float *z)
{
- double lon,lat,xp[3];
+ float lon,lat,xp[3];
int j,k;
/*
read in data for each layer
@@ -702,7 +714,7 @@
/*
read in lon lat
*/
- if(fscanf(in,"%lf %lf",&lon,&lat) != 2){
+ if(fscanf(in,"%f %f",&lon,&lat) != 2){
fprintf(stderr,"sh_read_spatial_data: error: lon lat format: pixel %i: read error\n",
(int)j);
exit(-1);
@@ -711,7 +723,7 @@
/*
read in lon lat z[i]
*/
- if(fscanf(in,"%lf %lf %f",&lon,&lat,z) != 3){
+ if(fscanf(in,"%f %f %f",&lon,&lat,z) != 3){
fprintf(stderr,"sh_read_spatial_data: error: lon lat z format: pixel %i: read error\n",
(int)j);
exit(-1);
@@ -776,7 +788,7 @@
hc_boolean verbose)
{
int j,os,inc;
- double xp[3];
+ float xp[3];
if(out_mode) /* make room for storing x,y,z */
hc_svecrealloc(x,exp->npoints*(2+((use_3d)?(1):(0))),"sh_compute_spatial_basis");
inc = (use_3d)?(3):(2);
@@ -882,7 +894,7 @@
*/
void sh_compute_spectral(float *data, int ivec,
- hc_boolean save_plm,double **plm,
+ hc_boolean save_plm,float **plm,
struct sh_lms *exp, hc_boolean verbose)
{
if(save_plm){
@@ -979,7 +991,7 @@
*/
void sh_compute_spatial(struct sh_lms *exp, int ivec,
- hc_boolean save_plm,double **plm,
+ hc_boolean save_plm,SH_RICK_PREC **plm,
float *data, hc_boolean verbose)
{
if((!exp[0].spectral_init)||(ivec && !exp[1].spectral_init)){
@@ -1001,7 +1013,7 @@
#ifdef HC_USE_HEALPIX
case SH_HEALPIX:
if(ivec){
- HC_ERROR("sh_compute_spectral","healpix: ivec==1 not implemented yet");
+ HC_ERROR("sh_compute_spatial","healpix: ivec==1 not implemented yet");
}else{
/*
scalar
@@ -1046,10 +1058,72 @@
sh_exp_type_error("sh_compute_spatial",exp);
break;
}
+}
+/*
-
+compute a spatial expansion on an irregular grid given in theta and
+phi arrays of npoints length
+
+cos(theta) and phi are cos(colatitude) and longitude in radians
+
+*/
+void sh_compute_spatial_irreg(struct sh_lms *exp, int ivec,
+ hc_boolean save_plm,SH_RICK_PREC **plm,
+ float *theta, int ntheta, float *phi,int nphi,
+ float *data, hc_boolean verbose,
+ hc_boolean save_sincos_fac)
+{
+ int npoints;
+ npoints = nphi * ntheta;
+ if((!exp[0].spectral_init)||(ivec && !exp[1].spectral_init)){
+ fprintf(stderr,"sh_compute_spatial_irreg: coefficients set not initialized, ivec: %i\n",
+ ivec);
+ exit(-1);
+ }
+ /*
+ make sure that the Legendre factors are computed
+ */
+ if(save_plm){
+ /*
+ this routine will only compute the plm once and performs
+ some sanity checks
+ */
+ sh_compute_plm_irreg(exp,ivec,plm,verbose,theta,ntheta);
+ }
+ switch(exp->type){
+#ifdef HC_USE_HEALPIX
+ case SH_HEALPIX:
+ HC_ERROR("sh_compute_spatial_irreg","healpix: not implemented yet");
+#endif
+ case SH_RICK:
+#ifdef NO_RICK_FORTRAN
+ if(save_plm)
+ rick_shc2d_pre_irreg(exp[0].alm,exp[1].alm,exp[0].lmax,
+ *plm,(*plm+exp->n_plm),
+ ivec,data,(data+npoints),
+ &exp->rick,theta,ntheta,phi,nphi,
+ save_sincos_fac);
+ else
+ rick_shc2d_irreg(exp[0].alm,exp[1].alm,exp[0].lmax,ivec,
+ data,(data+npoints),&exp->rick,
+ theta,ntheta,phi,nphi,save_sincos_fac);
+#else
+ HC_ERROR("sh_compute_spatial_irreg","Rick fortran not implemented");
+#endif
+ break;
+#ifdef HC_USE_SPHEREPACK
+ case SH_SPHEREPACK_GAUSS:
+ case SH_SPHEREPACK_EVEN:
+ HC_ERROR("sh_compute_spatial_irreg","Spherepack not implemented");
+ break;
+#endif
+ default:
+ sh_exp_type_error("sh_compute_spatial",exp);
+ break;
+ }
}
+
/*
print an error message and exit if a spherical harmonics type is not
@@ -1067,7 +1141,7 @@
print the Plm factors
*/
-void sh_print_plm(double *plm, int n_plm, int ivec, int type,
+void sh_print_plm(SH_RICK_PREC *plm, int n_plm, int ivec, int type,
FILE *out)
{
int i,j,jlim;
@@ -1120,7 +1194,7 @@
float z, FILE *out)
{
int j,k;
- double xp[3],lon,lat;
+ float xp[3],lon,lat;
for(j=0;j < exp[0].npoints;j++){
/*
get coordinates
@@ -1177,8 +1251,64 @@
fprintf(out,"\n");
} /* end points in lateral space loop */
}
+
/*
+irregular version
+
+ */
+void sh_print_irreg_spatial_data_to_file(struct sh_lms *exp, int shps,
+ float *data, hc_boolean use_3d,
+ float z, float *theta,int ntheta,
+ float *phi,int nphi,
+ FILE *out)
+{
+ int i,j,k,l,npoints;
+ float lon,lat;
+ npoints = nphi * ntheta;
+ /*
+ get coordinates
+ */
+ switch(exp->type){
+#ifdef HC_USE_HEALPIX
+
+ case SH_HEALPIX:
+ HC_ERROR("sh_print_irreg_spatial_data_to_file","healpix not implemented");
+ break;
+#endif
+ case SH_RICK:
+ for(i=l=00;i<ntheta;i++){
+ lat = THETA2LAT(theta[i]);
+ for(j=0;j<nphi;j++,l++){
+ lon = PHI2LON(phi[j]);
+ /* print coordinates */
+ if(!use_3d){
+ /* print lon lat */
+ fprintf(out,"%11g %11g\t",lon,lat);
+ }else{
+ /* print lon lat z[i] */
+ fprintf(out,"%11g %11g %11g\t",lon,lat,z);
+ }
+ for(k=0;k<shps;k++) /* loop through all scalars */
+ fprintf(out,"%11g ",data[l+npoints*k]);
+ fprintf(out,"\n");
+ }
+ }
+
+ break;
+#ifdef HC_USE_SPHEREPACK
+ case SH_SPHEREPACK_GAUSS:
+ case SH_SPHEREPACK_EVEN:
+ HC_ERROR("sh_print_irreg_spatial_data_to_file","spherepack not implemented");
+ break;
+#endif
+ default:
+ sh_exp_type_error("sh_print_spatial_data",exp);
+ break;
+ } /* end type branch */
+}
+/*
+
compute the associated Legendre functions for all (l,m) at
all latidutinal lcoations once and only once
@@ -1191,12 +1321,10 @@
plm: will be re-allocated, has to be passed at least as NULL
*/
-void sh_compute_plm(struct sh_lms *exp,int ivec,double **plm,
+void sh_compute_plm(struct sh_lms *exp,int ivec,SH_RICK_PREC **plm,
hc_boolean verbose)
{
- static hc_boolean plm_computed = FALSE;
- static int old_lmax,old_ivec,old_tnplm;
- if(!plm_computed){
+ if(!exp->plm_computed){
if((!exp->lmax)||(!exp->n_plm)||(!exp->tn_plm)){
fprintf(stderr,"sh_compute_plm: error, expansion not initialized?\n");
fprintf(stderr,"sh_compute_plm: lmax: %i n_plm: %i tn_plm: %i\n",
@@ -1206,7 +1334,7 @@
/*
allocate
*/
- hc_dvecrealloc(plm,exp->tn_plm,"sh_compute_plm");
+ hc_svecrealloc(plm,exp->tn_plm,"sh_compute_plm");
/*
compute the Legendre polynomials
*/
@@ -1241,36 +1369,121 @@
sh_exp_type_error("compute_plm",exp);
break;
}
- plm_computed = TRUE;
- old_lmax = exp->lmax;
- old_ivec = ivec;
- old_tnplm = exp->tn_plm;
+ exp->plm_computed = TRUE;
+ exp->old_lmax = exp->lmax;
+ exp->old_ivec = ivec;
+ exp->old_tnplm = exp->tn_plm;
}else{
/*
simple checks
first lmax
*/
- if(old_lmax != exp->lmax){
+ if(exp->old_lmax != exp->lmax){
fprintf(stderr,"sh_compute_plm: error: lmax initially %i, now %i\n",
- old_lmax,exp->lmax);
+ exp->old_lmax,exp->lmax);
exit(-1);
}
/* check if ivec was initialized if ever used */
- if(ivec > old_ivec){
+ if(ivec > exp->old_ivec){
fprintf(stderr,"sh_compute_plm: error: plm are to be saved but routine was initialized\n");
fprintf(stderr,"sh_compute_plm: error: with ivec: %i and now we want vectors, ivec: %i\n",
- old_ivec,ivec);
+ exp->old_ivec,ivec);
exit(-1);
}
- if(exp->tn_plm != old_tnplm){
+ if(exp->tn_plm != exp->old_tnplm){
fprintf(stderr,"sh_compute_plm: error: tn_plm initially %i, now %i\n",
- old_tnplm,exp->tn_plm);
+ exp->old_tnplm,exp->tn_plm);
exit(-1);
}
}
}
+
/*
+
+compute the associated Legendre functions for all (l,m) at all
+latidutinal lcoations once and only once for all irregularly placed latitudes in theya
+
+input:
+exp: holds the expansion parameters
+ivec_global: if 1, will construct vector arrays, else only for scalar
+
+output:
+
+plm: will be re-allocated, has to be passed at least as NULL
+
+*/
+void sh_compute_plm_irreg(struct sh_lms *exp,int ivec,SH_RICK_PREC **plm,
+ hc_boolean verbose, float *theta, int npoints)
+{
+ /* */
+ exp->tn_plm_irr = (1+ivec) * exp->lmsmall2 * npoints;
+
+ if((!exp->plm_computed_irr)||(exp->tn_plm_irr != exp->old_tnplm_irr)){
+ if((!exp->lmax)||(!exp->n_plm)||(!exp->tn_plm)){
+ fprintf(stderr,"sh_compute_plm_irreg: error, expansion not initialized?\n");
+ fprintf(stderr,"sh_compute_plm_irreg: lmax: %i n_plm: %i tn_plm: %i\n",
+ exp->lmax,exp->n_plm,exp->tn_plm);
+ exit(-1);
+ }
+ /*
+ allocate
+ */
+ exp->old_tnplm_irr = exp->tn_plm_irr;
+ hc_svecrealloc(plm,exp->old_tnplm_irr,"sh_compute_plm_irreg");
+ /*
+ compute the Legendre polynomials
+ */
+ switch(exp->type){
+#ifdef HC_USE_HEALPIX
+ HC_ERROR("compute_plm_irreg","healpix not implemented");
+ break;
+#endif
+ case SH_RICK:
+ if(verbose)
+ fprintf(stderr,"sh_compute_plm_irreg: Rick: computing all Plm for lmax %i and %i points\n",
+ exp->lmax,npoints);
+#ifdef NO_RICK_FORTRAN
+ rick_compute_allplm_irreg(exp->lmax,ivec,*plm,(*plm+exp->old_tnplm_irr),&exp->rick,
+ theta,npoints);
+#else
+ HC_ERROR("compute_plm_irreg","rick fortran not implemented");
+#endif
+ break;
+#ifdef HC_USE_SPHEREPACK
+ case SH_SPHEREPACK_GAUSS:
+ case SH_SPHEREPACK_EVEN:
+ HC_ERROR("compute_plm_irreg","spherepack implemented");
+ break;
+#endif
+ default:
+ sh_exp_type_error("compute_plm",exp);
+ break;
+ }
+ exp->plm_computed_irr = TRUE;
+ exp->old_lmax_irr = exp->lmax;
+ exp->old_ivec_irr = ivec;
+ }else{
+ /*
+ simple checks
+
+ first lmax
+ */
+ if(exp->old_lmax_irr != exp->lmax){
+ fprintf(stderr,"sh_compute_plm_irreg: error: lmax initially %i, now %i\n",
+ exp->old_lmax_irr,exp->lmax);
+ exit(-1);
+ }
+ /* check if ivec was initialized if ever used */
+ if(ivec > exp->old_ivec_irr){
+ fprintf(stderr,"sh_compute_plm_irreg: error: plm are to be saved but routine was initialized\n");
+ fprintf(stderr,"sh_compute_plm_irreg: error: with ivec: %i and now we want vectors, ivec: %i\n",
+ exp->old_ivec_irr,ivec);
+ exit(-1);
+ }
+ }
+}
+/*
returns the l,m A (use_b=0) or B (use_b=1) or A and B (use_b=2)
coefficient(s) of an expansion. if phys_normalization is set,
Modified: mc/1D/hc/trunk/sh_model.c
===================================================================
--- mc/1D/hc/trunk/sh_model.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/sh_model.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -72,8 +72,9 @@
up total number of points
in spatial domain
*/
+ /* use irregular grid */
sh_init_expansion((model->exp+i),lmax,type,model->ivec,
- verbose);
+ verbose,FALSE);
model->tnpoints += model->exp[i].npoints;
}
/* logic flag for spatial data */
Modified: mc/1D/hc/trunk/sh_power.c
===================================================================
--- mc/1D/hc/trunk/sh_power.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/sh_power.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -54,7 +54,7 @@
/*
input and init
*/
- sh_allocate_and_init(&exp,shps,lmax, type,ivec,verbose);
+ sh_allocate_and_init(&exp,shps,lmax, type,ivec,verbose,FALSE);
sh_read_coefficients_from_file(exp,shps,-1,stdin,binary,fac,verbose);
/* get space */
hc_svecrealloc(&power,exp->lmaxp1 * shps,"sh_power");
Modified: mc/1D/hc/trunk/sh_rick.h
===================================================================
--- mc/1D/hc/trunk/sh_rick.h 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/sh_rick.h 2007-07-12 19:08:51 UTC (rev 7650)
@@ -22,7 +22,7 @@
// how to convert from rick's spherical harmonic convention to dahlen and tromp?
// C_DT(l,m) = SH_RICK_FACTOR(l, m) * C_RICK(l,m)
//
-#define SH_RICK_FACTOR(l, m) (pow(-1.0,(double)(m))*SH_RICK_TWO_SQRT_PI)
+#define SH_RICK_FACTOR(l, m) (pow(-1.0,(SH_RICK_PREC)(m))*SH_RICK_TWO_SQRT_PI)
//
@@ -52,15 +52,18 @@
// Gauss points: cos(theta), weights, and actual theta
SH_RICK_PREC *gauss_z, *gauss_w, *gauss_theta;
//
+ //
+ SH_RICK_PREC *cfac,*sfac;
+
SH_RICK_PREC *lfac, *ilfac;
// those are for Legendre polynomials (fac1 fac2 only for ivec=1)
// make those double precision
//
- double *plm_f1,*plm_f2,*plm_fac1,*plm_fac2,*plm_srt;
+ SH_RICK_PREC *plm_f1,*plm_f2,*plm_fac1,*plm_fac2,*plm_srt;
// this is for vector harmonics, only for ivec=1
SH_RICK_PREC *sin_theta,*ell_factor;
// spacing in longitudes
- double dphi;
+ SH_RICK_PREC dphi;
// int (bounds and such)
int nlat,nlon,lmsize,lmsize2,nlonm1;
// logic flags
@@ -69,6 +72,10 @@
// init
my_boolean was_called;
+
+ /* more book keeping stuff */
+ int old_lmax,old_ivec,old_npoints,old_nplm,old_tnplm;
+
};
/*
@@ -104,19 +111,19 @@
extern void rick_f90_shd2c(float *,float *,int *,int *,
SH_RICK_PREC *,SH_RICK_PREC *);
-extern void rick_f90_shd2c_pre(float *, float *,int *,double *,
- double *, int *,SH_RICK_PREC *,
+extern void rick_f90_shd2c_pre(float *, float *,int *,SH_RICK_PREC *,
+ SH_RICK_PREC *, int *,SH_RICK_PREC *,
SH_RICK_PREC *);
extern void rick_f90_shc2d(SH_RICK_PREC *,SH_RICK_PREC *,int *,int *,
float *, float *);
extern void rick_f90_shc2d_pre(SH_RICK_PREC *,SH_RICK_PREC *,
- int *,double *,
- double *,int *,float *,float *);
+ int *,SH_RICK_PREC *,
+ SH_RICK_PREC *,int *,float *,float *);
extern void rick_f90_init(int *,int *,int *,int *, int *);
-extern void rick_f90_pix2ang(int *, int *, double *, double *);
+extern void rick_f90_pix2ang(int *, int *, SH_RICK_PREC *, SH_RICK_PREC *);
extern void rick_f90_free_module(int *);
extern void rick_f90_index(int *,int *,int *,int *);
-extern void rick_f90_compute_allplm(int *,int *,double *,double *);
+extern void rick_f90_compute_allplm(int *,int *,SH_RICK_PREC *,SH_RICK_PREC *);
/* FFT stuff */
extern void rick_f90_cs2ab( SH_RICK_PREC *, int *);
Modified: mc/1D/hc/trunk/sh_shana.h
===================================================================
--- mc/1D/hc/trunk/sh_shana.h 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/sh_shana.h 2007-07-12 19:08:51 UTC (rev 7650)
@@ -20,4 +20,6 @@
vector_sh_fac_init;
// init
hc_boolean was_called;
+
+ int old_ivec,old_lmax,old_npoints,old_tnplm;
};
Modified: mc/1D/hc/trunk/sh_syn.c
===================================================================
--- mc/1D/hc/trunk/sh_syn.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/sh_syn.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -16,11 +16,18 @@
int main(int argc, char **argv)
{
- int type,lmax,shps,ilayer,nset,ivec,i;
+ int type,lmax,shps,ilayer,nset,ivec,i,j,npoints,nphi,ntheta;
+ /*
+ switches
+ */
hc_boolean verbose = TRUE, short_format = FALSE ,short_format_ivec = FALSE ,binary = FALSE;
- float *data;
+ int regular_basis = 0;
+ /* */
+ float *data,*theta,*phi,x,y;
+ /* spacing for irregular output */
+ static float dphi = 2.0;
HC_PREC fac[3] = {1.,1.,1.},zlabel;
- double *dbl_dummy;
+ SH_RICK_PREC *dummy;
struct sh_lms *exp;
if(argc > 1){
if((strcmp(argv[1],"-h")==0)||(strcmp(argv[1],"--help")==0)||(strcmp(argv[1],"-help")==0))
@@ -32,17 +39,21 @@
}
}
if(argc > 2){
- sscanf(argv[1],"%i",&i);
+ sscanf(argv[2],"%i",&i);
if(i)
short_format_ivec = TRUE;
}
- if((argc > 3)||(argc<0)){
- fprintf(stderr,"usage: %s [short_format, %i] [short_ivec,%i]\n",
- argv[0],short_format,short_format_ivec);
+ if(argc > 3)
+ sscanf(argv[3],"%i",®ular_basis);
+
+ if((argc > 4)||(argc<0)){
+ fprintf(stderr,"usage: %s [short_format, %i] [short_ivec,%i] [regular_basis,%i]\n",
+ argv[0],short_format,short_format_ivec,regular_basis);
exit(-1);
}
- fprintf(stderr,"%s: waiting to read spherical harmonic coefficients from stdin\n",
- argv[0]);
+ if(verbose)
+ fprintf(stderr,"%s: waiting to read spherical harmonic coefficients from stdin, use -h for help\n",
+ argv[0]);
while(sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer,&nset,
&zlabel,&ivec,stdin,short_format,
binary,verbose)){
@@ -50,17 +61,58 @@
ivec = 1;
shps = 2;
}
- fprintf(stderr,"%s: converting lmax %i ivec: %i at z: %g\n",
- argv[0],lmax,ivec,zlabel);
+ if(verbose)
+ fprintf(stderr,"%s: converting lmax %i ivec: %i at z: %g\n",
+ argv[0],lmax,ivec,zlabel);
/* input and init */
- sh_allocate_and_init(&exp,shps,lmax,type,ivec,verbose);
+ sh_allocate_and_init(&exp,shps,lmax,type,ivec,verbose,regular_basis);
sh_read_coefficients_from_file(exp,shps,-1,stdin,binary,fac,verbose);
- /* expansion */
- hc_svecalloc(&data,exp[0].npoints * shps,"sh_shsyn");
- sh_compute_spatial(exp,ivec,FALSE,&dbl_dummy,data,verbose);
- /* output */
- sh_print_spatial_data_to_file(exp,shps,data,(nset>1)?(TRUE):(FALSE),zlabel,stdout);
+ if(regular_basis){
+ /*
+ irregular basis output
+ */
+ switch(regular_basis){
+ case 1:
+ if(verbose)
+ fprintf(stderr,"sh_syn: using regular spaced grid with %g deg spacing\n",dphi);
+ /* */
+ dphi = DEG2RAD(dphi);
+ nphi = (TWOPI-dphi)/dphi;
+ ntheta = nphi/2+1;
+ npoints = nphi * ntheta;
+ /* */
+ hc_svecalloc(&phi,nphi,"sh_shsyn");
+ hc_svecalloc(&theta,ntheta,"sh_shsyn");
+ for(x=0,i=0;i < nphi;i++,x+=dphi)
+ phi[i] = x;
+ for(y=dphi/2,j=0;j<ntheta;y+=dphi,j++)
+ theta[j] = y;
+ break;
+ default:
+ fprintf(stderr,"sh_syn: irregular basis code %i undefined\n",regular_basis);
+ exit(-1);
+ break;
+ }
+
+ hc_svecalloc(&data,npoints * shps,"sh_shsyn data");
+ sh_compute_spatial_irreg(exp,ivec,FALSE,&dummy,
+ theta,ntheta,phi,nphi,data,verbose,FALSE);
+ /* output */
+ sh_print_irreg_spatial_data_to_file(exp,shps,data,(nset>1)?(TRUE):(FALSE),zlabel,
+ theta,ntheta,phi,nphi,stdout);
+
+
+
+ }else{ /* use the built in spatial basis (Gaussian) */
+ /* expansion */
+ hc_svecalloc(&data,exp[0].npoints * shps,"sh_syn");
+ sh_compute_spatial(exp,ivec,FALSE,&dummy,data,verbose);
+ /* output */
+ sh_print_spatial_data_to_file(exp,shps,data,(nset>1)?(TRUE):(FALSE),zlabel,stdout);
+ }
+
+
free(exp);free(data);
}
Modified: mc/1D/hc/trunk/shana_sh.c
===================================================================
--- mc/1D/hc/trunk/shana_sh.c 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/shana_sh.c 2007-07-12 19:08:51 UTC (rev 7650)
@@ -498,7 +498,6 @@
// local
HC_CPREC xtemp;
int i,l;
- static int old_lmax,old_ivec,old_npoints,old_nplm,old_tnplm;
if(!shana->was_called){
if(lmax == 0){
@@ -525,7 +524,7 @@
// number of points in one layer
//
*npoints = shana->nlat * shana->nlon;
- old_npoints = *npoints;
+ shana->old_npoints = *npoints;
//
//
// for coordinate computations
@@ -541,8 +540,8 @@
// size of the Plm array
*nplm = shana->lmsize * shana->nlat;
*tnplm = *nplm * (1+ivec); // for all layers
- old_tnplm = *tnplm;
- old_nplm = *nplm;
+ shana->old_tnplm = *tnplm;
+ shana->old_nplm = *nplm;
//
// initialize the Gauss points, at which the latitudes are
// evaluated
@@ -596,23 +595,23 @@
/*
save initial call lmax and ivec settings
*/
- old_lmax = lmax;
- old_ivec = ivec;
+ shana->old_lmax = lmax;
+ shana->old_ivec = ivec;
/* end initial branch */
}else{
- if(lmax != old_lmax){
+ if(lmax != shana->old_lmax){
fprintf(stderr,"shana_init: error: was init with lmax %i, now: %i. (ivec: %i, now: %i)\n",
- old_lmax,lmax,old_ivec,ivec);
+ shana->old_lmax,lmax,shana->old_ivec,ivec);
exit(-1);
}
- if(ivec > old_ivec){
- fprintf(stderr,"shana_init: error: original ivec %i, now %i\n",old_ivec,ivec);
+ if(ivec > shana->old_ivec){
+ fprintf(stderr,"shana_init: error: original ivec %i, now %i\n",shana->old_ivec,ivec);
exit(-1);
}
- *npoints = old_npoints;
- *nplm = old_nplm;
- *tnplm = old_tnplm;
+ *npoints = shana->old_npoints;
+ *nplm = shana->old_nplm;
+ *tnplm = shana->old_tnplm;
}
} /* end shana init */
@@ -645,7 +644,7 @@
double plm,pm1,pm2,pmm,sintsq,fnum,fden;
//
int i,l,m,k,kstart,l2,mstop,lmaxm1;
- static int old_lmax,old_ivec;
+
if(!shana->initialized){
fprintf(stderr,"shana_plmbar1: error: module not initialized, call shana_init first\n");
exit(-1);
@@ -716,8 +715,8 @@
k++;
}
} /* end ivec==1 */
- old_lmax = lmax;
- old_ivec = ivec;
+ shana->old_lmax = lmax;
+ shana->old_ivec = ivec;
shana->computed_legendre = TRUE;
/*
end first call
@@ -727,13 +726,13 @@
subsequent call
*/
// test if lmax has changed
- if(lmax != old_lmax){
- fprintf(stderr,"shana_plmbar1: error: factors were computed for lmax %in",old_lmax);
+ if(lmax != shana->old_lmax){
+ fprintf(stderr,"shana_plmbar1: error: factors were computed for lmax %in",shana->old_lmax);
fprintf(stderr,"shana_plmbar1: error: now, lmax is %i\n",lmax);
exit(-1);
}
- if(ivec > old_ivec){
- fprintf(stderr,"shana_plmbar1: error: init with %i, now ivec %i\n",old_ivec,ivec);
+ if(ivec > shana->old_ivec){
+ fprintf(stderr,"shana_plmbar1: error: init with %i, now ivec %i\n",shana->old_ivec,ivec);
exit(-1);
}
}
Modified: mc/1D/hc/trunk/test_sh_routines
===================================================================
--- mc/1D/hc/trunk/test_sh_routines 2007-07-12 18:33:54 UTC (rev 7649)
+++ mc/1D/hc/trunk/test_sh_routines 2007-07-12 19:08:51 UTC (rev 7650)
@@ -3,9 +3,9 @@
# test spherical harmonic routines for scalar fields
#
makegrid=${1-0}
-vec=0 # 0: scalar test 1: velocity field test
+vec=1 # 0: scalar test 1: velocity field test
lmax=31
-type=1 # 0: Healpix 1: Rick
+type=0 # 1: Healpix 0: Rick
More information about the cig-commits
mailing list