[cig-commits] r7893 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Sun Aug 26 21:43:03 PDT 2007
Author: becker
Date: 2007-08-26 21:43:01 -0700 (Sun, 26 Aug 2007)
New Revision: 7893
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_grdtrack_util.c
mc/1D/hc/trunk/ggrd_grdtrack_util.h
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_extract_sh_layer.c
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/rick_sh_c.c
mc/1D/hc/trunk/sh_ana.c
mc/1D/hc/trunk/sh_exp.c
mc/1D/hc/trunk/sh_rick.h
mc/1D/hc/trunk/sh_syn.c
mc/1D/hc/trunk/vdepth.dat
Log:
- Fixed a few minor bugs.
- Added a new example script for running hc, and plotting flow solutions
using GMT: calc_vel_and_plot
Modified: mc/1D/hc/trunk/.gmtcommands
===================================================================
--- mc/1D/hc/trunk/.gmtcommands 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/.gmtcommands 2007-08-27 04:43:01 UTC (rev 7893)
@@ -1,3 +1,7 @@
# GMT common arguments shelf
+-B1/:v at -r@- [cm/yr]:
+-JH180/7
+-JX7/3.5
-R0/360/-90/90
EOF
+EOF
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/Makefile 2007-08-27 04:43:01 UTC (rev 7893)
@@ -9,7 +9,7 @@
#
#
LD = $(CC)
-CFLAGS = $(CFLAGS_DEBUG)
+#CFLAGS = $(CFLAGS_DEBUG)
#
# EDIT HERE FOR GMT VERSION
Modified: mc/1D/hc/trunk/README.TXT
===================================================================
--- mc/1D/hc/trunk/README.TXT 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/README.TXT 2007-08-27 04:43:01 UTC (rev 7893)
@@ -194,7 +194,11 @@
+ Also see the script calc_vel_and_plot for some suggestions on
+ how to convert the output.
+
+
USING THE SPHERICAL HARMONICS TOOLS AS STANDALONE
Modified: mc/1D/hc/trunk/expansion.par
===================================================================
--- mc/1D/hc/trunk/expansion.par 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/expansion.par 2007-08-27 04:43:01 UTC (rev 7893)
@@ -1 +1 @@
- 0.0000000000 2.0000000000 360.0000000000 -90.0000000000 2.0000000000 90.0000000000 181 91
+ 0.0000000000 1.0000000000 360.0000000000 -90.0000000000 1.0000000000 90.0000000000 361 181
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2007-08-27 04:43:01 UTC (rev 7893)
@@ -832,7 +832,7 @@
fprintf(stderr,"ggrd_read_time_intervals: error: already initialized\n");
return 1;
}
- hc_vecalloc(&thist->vtimes,3,"rti: 1");
+ ggrd_vecalloc(&thist->vtimes,3,"rti: 1");
if(read_thistory){
in = fopen(input_file,"r");
@@ -859,7 +859,7 @@
(*(thist->vtimes+thist->nvtimes3) + *(thist->vtimes + thist->nvtimes3+2))/2.0;
thist->nvtimes += 1;
thist->nvtimes3 += 3;
- hc_vecrealloc(&thist->vtimes,thist->nvtimes3+3,"rti: 2");
+ ggrd_vecrealloc(&thist->vtimes,thist->nvtimes3+3,"rti: 2");
}
thist->tmin = *(thist->vtimes+0);
thist->tmax = *(thist->vtimes+ (thist->nvtimes-1) * 3 +2);
@@ -1161,6 +1161,29 @@
return 0;
}
+
+
+/* general floating point vector allocation */
+void ggrd_vecalloc(double **x,int n,char *message)
+{
+ *x = (double *)malloc(sizeof(double)*(size_t)n);
+ if(! (*x)){
+ fprintf(stderr,"mem error: %s\n",message);
+ exit(-1);
+ }
+}
+
+/* general version */
+void ggrd_vecrealloc(double **x,int n,char *message)
+{
+ *x = (double *)realloc(*x,sizeof(double)*(size_t)n);
+ if(!(*x)){
+ fprintf(stderr,"mem error: %s\n",message);
+ exit(-1);
+ }
+}
+
+
/*
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.h
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.h 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.h 2007-08-27 04:43:01 UTC (rev 7893)
@@ -103,7 +103,6 @@
double *);
-/* hc related utility functions */
-void hc_vecalloc(double **,int,char *);
-void hc_vecrealloc(double **,int,char *);
+void ggrd_vecalloc(double **,int,char *);
+void ggrd_vecrealloc(double **,int,char *);
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/hc.h 2007-08-27 04:43:01 UTC (rev 7893)
@@ -11,6 +11,7 @@
#include <math.h>
#include <stdlib.h>
#include <string.h>
+#include <malloc.h>
#include <limits.h>
#include <malloc.h>
@@ -131,6 +132,10 @@
hc_boolean print_spatial; /* print the spatial solution */
hc_boolean compute_geoid; /* compute and print the geoid */
int solution_mode; /* velocity or stress */
+
+ hc_boolean read_short_dens_sh; /* short SH format for density
+ files? */
+
hc_boolean print_pt_sol; /* output of p[6] and t[2] vectors */
char visc_filename[HC_CHAR_LENGTH]; /* name of viscosity profile file */
char pvel_filename[HC_CHAR_LENGTH]; /* name of plate velocities file */
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2007-08-27 04:43:01 UTC (rev 7893)
@@ -12,6 +12,8 @@
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 *);
+void ggrd_vecalloc(double **, int, char *);
+void ggrd_vecrealloc(double **, int, char *);
float ggrd_gt_rms(float *, int);
float ggrd_gt_mean(float *, int);
/* ggrd_readgrds.c */
@@ -33,7 +35,7 @@
void hc_init_constants(struct hcs *, double, char *, unsigned short);
void hc_handle_command_line(int, char **, struct hc_parameters *);
void hc_assign_viscosity(struct hcs *, int, char [300], unsigned short);
-void hc_assign_density(struct hcs *, unsigned short, int, char *, int, unsigned short, unsigned short, unsigned short);
+void hc_assign_density(struct hcs *, unsigned short, int, char *, int, unsigned short, unsigned short, unsigned short, unsigned short);
void hc_init_phase_boundaries(struct hcs *, int, unsigned short);
void hc_assign_plate_velocities(struct hcs *, int, char *, unsigned short, int, unsigned short, unsigned short);
void hc_init_l_factors(struct hcs *, int);
Modified: mc/1D/hc/trunk/hc_extract_sh_layer.c
===================================================================
--- mc/1D/hc/trunk/hc_extract_sh_layer.c 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/hc_extract_sh_layer.c 2007-08-27 04:43:01 UTC (rev 7893)
@@ -53,10 +53,13 @@
fprintf(stderr,"\tif mode = 1, will print x_r \n");
fprintf(stderr,"\tif mode = 2, will print x_pol x_tor \n");
fprintf(stderr,"\tif mode = 3, will print x_r x_pol x_tor\n");
+ fprintf(stderr,"\tif mode = 4, will print the depth levels of all layers\n");
exit(-1);
break;
}
+ if(mode == 4)
+ ilayer = -2;
/*
read in solution
*/
@@ -92,11 +95,14 @@
/*
output
*/
- if(short_format && loop)
- fprintf(stdout,"%g\n",HC_Z_DEPTH(model->r[ilayer]));
- sh_print_parameters_to_file((sol+ilayer*3),shps,
- ilayer,nset,(HC_PREC)(HC_Z_DEPTH(model->r[ilayer])),
- stdout,short_format,FALSE,verbose);
+ if(mode != 4){
+ /* SH header */
+ if(short_format && loop)
+ fprintf(stdout,"%g\n",HC_Z_DEPTH(model->r[ilayer]));
+ sh_print_parameters_to_file((sol+ilayer*3),shps,
+ ilayer,nset,(HC_PREC)(HC_Z_DEPTH(model->r[ilayer])),
+ stdout,short_format,FALSE,verbose);
+ }
switch(mode){
case 1:
/* */
@@ -119,6 +125,9 @@
argv[0],ilayer,HC_Z_DEPTH(model->r[ilayer]));
sh_print_coefficients_to_file((sol+ilayer*3),shps,stdout,fac,FALSE,verbose);
break;
+ case 4:
+ fprintf(stdout,"%5i %11g\n",ilayer,HC_Z_DEPTH(model->r[ilayer]));
+ break;
default:
fprintf(stderr,"%s: error, mode %i undefined\n",argv[0],mode);
exit(-1);
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/hc_init.c 2007-08-27 04:43:01 UTC (rev 7893)
@@ -29,6 +29,10 @@
p->dens_anom_scale = 0.2; /* default density anomaly scaling to
go from PREM percent traveltime
anomalies to density anomalies */
+ p->read_short_dens_sh = FALSE; /* read the density anomaly file in
+ the short format of Becker &
+ Boschi (2002)? */
+
p->verbose = 0; /* debugging output? (0,1,2,3,4...) */
p->sol_binary_out = TRUE; /* binary or ASCII output of SH expansion */
p->solution_mode = HC_VEL; /* default: velocity */
@@ -121,7 +125,9 @@
HC_ERROR("hc_init","vel_bc_zero and free_slip doesn't make sense");
/* read in the densities */
hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,
- p->dens_filename,-1,FALSE,FALSE,p->verbose);
+ p->dens_filename,-1,FALSE,FALSE,
+ p->verbose,
+ p->read_short_dens_sh);
/* assign all zeroes up to the lmax of the density expansion */
hc_assign_plate_velocities(hc,HC_INIT_FROM_FILE,
p->pvel_filename,
@@ -138,14 +144,15 @@
dummy,FALSE,p->verbose);
hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,
p->dens_filename,hc->pvel[0].lmax,
- FALSE,FALSE,p->verbose);
+ FALSE,FALSE,p->verbose,
+ p->read_short_dens_sh);
}else{
if(p->verbose)
fprintf(stderr,"hc_init: initializing for free-slip\n");
/* read in the density fields */
hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,
p->dens_filename,-1,FALSE,FALSE,
- p->verbose);
+ p->verbose,p->read_short_dens_sh);
}
}
/*
@@ -250,7 +257,7 @@
{
int i;
for(i=1;i < argc;i++){
- if(strcmp(argv[i],"-h")==0 || strcmp(argv[i],"-?")==0){// help
+ if(strcmp(argv[i],"-help")==0 || strcmp(argv[i],"-h")==0 || strcmp(argv[i],"-?")==0){// help
/*
help page
*/
@@ -264,7 +271,11 @@
fprintf(stderr,"density anomaly options:\n");
fprintf(stderr,"-dens\tname\tuse name as a SH density anomaly model (%s)\n",
p->dens_filename);
- fprintf(stderr,"\t\tThis expects density anomalies in %% of PREM in DT format.\n");
+ fprintf(stderr,"\t\tThis expects density anomalies in %% of PREM in Dahlen & Tromp format.\n");
+
+ fprintf(stderr,"-dshs\t\tuse the short, Becker & Boschi (2002) format for the SH density model (%s)\n",
+ hc_name_boolean(p->read_short_dens_sh));
+
fprintf(stderr,"-ds\t\tadditional density anomaly scaling factor (%g)\n\n",
p->dens_anom_scale);
@@ -305,6 +316,8 @@
hc_toggle_boolean(&p->vel_bc_zero);
}else if(strcmp(argv[i],"-px")==0){ /* print spatial solution? */
hc_toggle_boolean(&p->print_spatial);
+ }else if(strcmp(argv[i],"-dshs")==0){ /* use short format? */
+ hc_toggle_boolean(&p->read_short_dens_sh);
}else if(strcmp(argv[i],"-pptsol")==0){ /* print
poloidal/toroidal
solution
@@ -477,12 +490,14 @@
char *filename,int nominal_lmax,
hc_boolean layer_structure_changed,
hc_boolean density_in_binary,
- hc_boolean verbose)
+ hc_boolean verbose,
+ hc_boolean use_short_format)
{
FILE *in;
int type,lmax,shps,ilayer,nset,ivec,i,j;
HC_PREC *dtop,*dbot,zlabel,dens_scale[1],rho0;
- hc_boolean reported = FALSE;
+ double dtmp;
+ hc_boolean reported = FALSE,read_on;
static HC_PREC local_dens_fac = .01; /* this factor will be
multiplied with the
hc->dens_fac factor to
@@ -492,7 +507,7 @@
for instance
*/
hc->compressible = compressible;
-
+ hc->inho = 0;
if(hc->dens_init) /* clear old expansions, if
already initialized */
sh_free_expansion(hc->dens_anom,hc->inho);
@@ -512,98 +527,128 @@
this assumes that we are reading in anomalies in percent
*/
+
in = hc_open(filename,"r","hc_assign_density");
if(verbose)
fprintf(stderr,"hc_assign_density: reading density anomalies in [%%] from %s, scaling by %g\n",
filename,hc->dens_scale[0]);
hc->inho = 0; /* counter for density layers */
/* get one density expansion */
- hc_get_blank_expansions(&hc->dens_anom,1,0,"hc_assign_density");
+ hc_get_blank_expansions(&hc->dens_anom,1,0,
+ "hc_assign_density");
/*
read all layers as spherical harmonics assuming real Dahlen &
Tromp (physical) normalization, short format
*/
- while(sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer, &nset,
- &zlabel,&ivec,in,FALSE,
- density_in_binary,
- verbose)){
- if((verbose)&&(!reported)){
- if(nominal_lmax > lmax)
- fprintf(stderr,"hc_assign_density: density lmax: %3i filling up to nominal lmax: %3i with zeroes\n",
- lmax,nominal_lmax);
- if(nominal_lmax != -1){
- fprintf(stderr,"hc_assign_density: density lmax: %3i limiting to lmax: %3i\n",
- lmax,nominal_lmax);
- }else{
- fprintf(stderr,"hc_assign_density: density lmax: %3i determines solution lmax\n",
- lmax);
+ if(use_short_format){
+ if(verbose)
+ fprintf(stderr,"hc_assign_density: using short format for density SH\n");
+ fscanf(in,"%i",&nset);
+ ilayer = -1;
+ }else{
+ if(verbose)
+ fprintf(stderr,"hc_assign_density: using default SH format for density\n");
+ }
+ read_on = TRUE;
+ while(read_on){
+ if(use_short_format){
+ /* short format I/O */
+ i = fscanf(in,"%lf",&dtmp);zlabel = (HC_PREC)dtmp;
+ i += fscanf(in,"%i",&lmax);
+ read_on = (i == 2)?(TRUE):(FALSE);
+ ivec = 0;shps = 1;type = HC_DEFAULT_INTERNAL_FORMAT;
+ ilayer++;
+ }else{
+ read_on = sh_read_parameters_from_file(&type,&lmax,&shps,
+ &ilayer, &nset,
+ &zlabel,&ivec,in,
+ FALSE,density_in_binary,
+ verbose);
+ }
+ if(read_on){
+ if((verbose)&&(!reported)){
+ if(nominal_lmax > lmax)
+ fprintf(stderr,"hc_assign_density: density lmax: %3i filling up to nominal lmax: %3i with zeroes\n",
+ lmax,nominal_lmax);
+ if(nominal_lmax != -1){
+ fprintf(stderr,"hc_assign_density: density lmax: %3i limiting to lmax: %3i\n",
+ lmax,nominal_lmax);
+ }else{
+ fprintf(stderr,"hc_assign_density: density lmax: %3i determines solution lmax\n",
+ lmax);
+ }
+ reported = TRUE;
}
- reported = TRUE;
- }
- /*
- do tests
- */
- if((shps != 1)||(ivec))
- HC_ERROR("hc_assign_density","vector field read in but only scalar expansion expected");
- /* test and assign depth levels */
- hc->rden=(HC_PREC *)
- realloc(hc->rden,(1+hc->inho)*sizeof(HC_PREC));
- if(!hc->rden)
- HC_MEMERROR("hc_assign_density: rden");
- /*
- assign depth, this assumes that we are reading in depths [km]
- */
- hc->rden[hc->inho] = HC_ND_RADIUS((double)zlabel);
- /*
-
- get reference density at this level
-
- */
- prem_get_rho(&rho0,hc->rden[hc->inho],hc->prem);
- rho0 /= 1000.0;
- /*
- general density (add additional depth dependence here)
- */
- dens_scale[0] = hc->dens_scale[0] * local_dens_fac * rho0;
- if(verbose >= 2)
- fprintf(stderr,"hc_assign_density: r: %11g anom scales: %11g x %11g x %11g = %11g\n",
- hc->rden[hc->inho],hc->dens_scale[0],
- local_dens_fac,rho0,dens_scale[0]);
- if(hc->inho){
/*
- check by comparison with previous expansion
+ do tests
*/
- if(nominal_lmax == -1)
- if(lmax != hc->dens_anom[0].lmax)
- HC_ERROR("hc_assign_density","lmax changed in file");
- if(hc->rden[hc->inho] <= hc->rden[hc->inho-1])
- HC_ERROR("hc_assign_density","depth should decrease, radius increase (give z[km])");
- }
- /*
- make room for new expansion
- */
- hc_get_blank_expansions(&hc->dens_anom,hc->inho+1,hc->inho,
- "hc_assign_density");
- /*
- initialize expansion on irregular grid
- */
- sh_init_expansion((hc->dens_anom+hc->inho),
- (nominal_lmax > lmax) ? (nominal_lmax):(lmax),
- hc->sh_type,1,verbose,FALSE);
- /*
-
- read parameters and scale (put possible depth dependence of
- scaling here)
-
- will assume input is in physical convention
-
- */
- sh_read_coefficients_from_file((hc->dens_anom+hc->inho),1,lmax,
- in,density_in_binary,dens_scale,
- verbose);
- hc->inho++;
- }
+ if((shps != 1)||(ivec))
+ HC_ERROR("hc_assign_density","vector field read in but only scalar expansion expected");
+ /* test and assign depth levels */
+ hc->rden=(HC_PREC *)
+ realloc(hc->rden,(1+hc->inho)*sizeof(HC_PREC));
+ if(!hc->rden)
+ HC_MEMERROR("hc_assign_density: rden");
+ /*
+ assign depth, this assumes that we are reading in depths [km]
+ */
+ hc->rden[hc->inho] = HC_ND_RADIUS((double)zlabel);
+ /*
+
+ get reference density at this level
+
+ */
+ prem_get_rho(&rho0,hc->rden[hc->inho],hc->prem);
+ rho0 /= 1000.0;
+ /*
+ general density (add additional depth dependence here)
+ */
+ dens_scale[0] = hc->dens_scale[0] * local_dens_fac * rho0;
+ if(verbose >= 2){
+
+ fprintf(stderr,"hc_assign_density: r: %11g anom scales: %11g x %11g x %11g = %11g\t%5i out of %i, z: %11g\n",
+ hc->rden[hc->inho],hc->dens_scale[0],
+ local_dens_fac,rho0,dens_scale[0],hc->inho+1,nset,zlabel);
+ }
+ if(hc->inho){
+ /*
+ check by comparison with previous expansion
+ */
+ if(nominal_lmax == -1)
+ if(lmax != hc->dens_anom[0].lmax)
+ HC_ERROR("hc_assign_density","lmax changed in file");
+ if(hc->rden[hc->inho] <= hc->rden[hc->inho-1]){
+ fprintf(stderr,"hc_assign_density: %i %g %g\n",hc->inho,hc->rden[hc->inho] <= hc->rden[hc->inho-1]);
+ HC_ERROR("hc_assign_density","depth should decrease, radius increase (give z[km])");
+ }
+ }
+ /*
+ make room for new expansion
+ */
+ hc_get_blank_expansions(&hc->dens_anom,hc->inho+1,hc->inho,
+ "hc_assign_density");
+ /*
+ initialize expansion on irregular grid
+ */
+ sh_init_expansion((hc->dens_anom+hc->inho),
+ (nominal_lmax > lmax) ? (nominal_lmax):(lmax),
+ hc->sh_type,1,verbose,FALSE);
+ /*
+
+ read parameters and scale (put possible depth dependence of
+ scaling here)
+
+ will assume input is in physical convention
+
+ */
+ sh_read_coefficients_from_file((hc->dens_anom+hc->inho),1,
+ lmax,in,density_in_binary,
+ dens_scale,verbose);
+ hc->inho++;
+ } /* end actualy read on */
+ } /* end while */
+
if(hc->inho != nset)
HC_ERROR("hc_assign_density","file mode: mismatch in number of layers");
fclose(in);
@@ -884,8 +929,6 @@
{
free((*hc)->visc);
free((*hc)->rvisc);
- free((*hc)->visc);
-
free((*hc)->qwrite);
sh_free_expansion((*hc)->dens_anom,1);
Modified: mc/1D/hc/trunk/rick_sh_c.c
===================================================================
--- mc/1D/hc/trunk/rick_sh_c.c 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/rick_sh_c.c 2007-08-27 04:43:01 UTC (rev 7893)
@@ -128,12 +128,13 @@
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)
+ float *phi,int nphi,
+ hc_boolean save_sincos_fac)
{
SH_RICK_PREC *plm,*dplm;
if(!rick->initialized){
@@ -316,7 +317,8 @@
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)
+ float *phi,int nphi,
+ my_boolean save_sincos_fac)
{
/* //
// Legendre functions are precomputed
@@ -329,7 +331,7 @@
exit(-1);
}
- lmaxp1 = lmax + 1; // this is nlat
+ lmaxp1 = lmax + 1;
if(ivec){
if(!rick->vector_sh_fac_init){
@@ -345,14 +347,19 @@
/*
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((!save_sincos_fac)||(!rick->sin_cos_saved)){
+
+ 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(save_sincos_fac)
+ rick->sin_cos_saved = TRUE;
}
if(ivec == 0){
/*
@@ -368,8 +375,7 @@
loc_plmb[ios1] = cslm[k2+1] * plm[ios1];
}
}
- for (idata=i=ios2=0;i < ntheta; i++,ios2 += rick->lmsize) { /* theta
- loop */
+ 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*/
@@ -397,17 +403,22 @@
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;
+ l = 0;m = 0;
/* start at l = 1 */
for (k=1,k2=2; k < rick->lmsize; k++,k2+=2) {
- /* loop through l,m */
+ /* loop through l,m */
m++;
if (m > l) {
m=0;l++;
}
lm1 = l - 1;
+ // fprintf(stderr,"%5i %5i\t %11g %11g\tPA: %11g PB: %11g\tTA: %11g TB:%11g\n",
+ // l,m,rick->ell_factor[lm1],plm[ios2+k],
+ // SH_RICK_FACTOR(l, m) * cslm[k2] ,SH_RICK_FACTOR(l, m) *cslm[k2+1],
+ // SH_RICK_FACTOR(l, m) * dslm[k2], SH_RICK_FACTOR(l, m) *dslm[k2+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 */
@@ -418,21 +429,23 @@
u_theta
- */
- rdatax[i] += /* cos term */
+ */
+ rdatax[idata] += /* cos term */
(cslm[k2] * dpdt + dslm[k2+1] * dpdp) * rick->cfac[ios3+m];
- rdatax[i] += /* sin term */
+ rdatax[idata] += /* sin term */
(cslm[k2+1] * dpdt - dslm[k2] * dpdp) * rick->sfac[ios3+m];
/*
u_phi
*/
- rdatay[i] += // cos term
+ rdatay[idata] += // cos term
(cslm[k2+1] * dpdp - dslm[k2] * dpdt) * rick->cfac[ios3+m];
- rdatay[i] += // sin term
+ rdatay[idata] += // sin term
(- cslm[k2] * dpdp - dslm[k2+1] * dpdt) * rick->sfac[ios3+m];
}
- }
- }
+ //fprintf(stderr,"%11g %11g %11g %11g\n",theta[j],phi[i],rdatax[i],rdatay[i]);
+
+ } /* end phi loop */
+ } /* end theta loop */
}
if(!save_sincos_fac){
@@ -759,6 +772,10 @@
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");
+
+
+ rick->sin_cos_saved = FALSE;
+
if(ivec){
//
// additional arrays for vector spherical harmonics
Modified: mc/1D/hc/trunk/sh_ana.c
===================================================================
--- mc/1D/hc/trunk/sh_ana.c 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/sh_ana.c 2007-08-27 04:43:01 UTC (rev 7893)
@@ -38,7 +38,7 @@
int main(int argc, char **argv)
{
- int type = SH_RICK,lmax,shps, nset=1,ivec=0,ilayer=0;
+ int type = SH_RICK,lmax,shps, nset=1,ivec=0,ilayer=0,i;
struct sh_lms *exp;
hc_boolean verbose = TRUE, use_3d = FALSE, short_format = FALSE,
binary = FALSE, print_spatial_base = FALSE;
@@ -64,16 +64,22 @@
sscanf(argv[2],"%i",&ivec);
if(argc > 3)
sscanf(argv[3],"%i",&type);
+ if(argc > 4){
+ sscanf(argv[4],"%i",&i);
+ short_format = (hc_boolean)i;
+ }
if((argc > 4)||(argc<=1)){
- fprintf(stderr,"usage: %s l_max [ivec, %i] [type, %i]\n",
- argv[0],ivec,type);
- fprintf(stderr," l_max: max order of expansion. if negative, will print out the spatial\n");
- fprintf(stderr," locations needed on input\n");
- fprintf(stderr," for Rick, lmax needs to be 2**n-1\n\n");
- fprintf(stderr," ivec: 0: expand scalar field (input: lon lat scalar)\n");
- fprintf(stderr," 1: expand vector field (input: lon lat v_t v_p\n\n");
- fprintf(stderr," type: %i: use Healpix's routines\n",SH_HEALPIX);
- fprintf(stderr," %i: use Rick's routines\n",SH_RICK);
+ fprintf(stderr,"usage: %s l_max [ivec, %i] [type, %i] [short_format, %i\n",
+ argv[0],ivec,type,short_format);
+ fprintf(stderr," l_max: max order of expansion. if negative, will print out the spatial\n");
+ fprintf(stderr," locations needed on input\n");
+ fprintf(stderr," for Rick SH format, lmax needs to be 2**n-1\n\n");
+ fprintf(stderr," ivec: 0: expand scalar field (input: lon lat scalar)\n");
+ fprintf(stderr," 1: expand vector field (input: lon lat v_t v_p\n\n");
+ fprintf(stderr," type %i: use Rick's routines internally (output format DT convection)\n",SH_RICK);
+ fprintf(stderr," %i: use Healpix's routines internally (output format DT convection)\n\n",SH_HEALPIX);
+ fprintf(stderr,"short_format: 0: use long header files format as in HC\n");
+ fprintf(stderr," 1: use short header file format (only lmax)\n\n");
exit(-1);
}
if(print_spatial_base)
Modified: mc/1D/hc/trunk/sh_exp.c
===================================================================
--- mc/1D/hc/trunk/sh_exp.c 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/sh_exp.c 2007-08-27 04:43:01 UTC (rev 7893)
@@ -363,11 +363,15 @@
*/
-hc_boolean sh_read_parameters_from_file(int *type, int *lmax, int *shps,
+hc_boolean sh_read_parameters_from_file(int *type, int *lmax,
+ int *shps,
int *ilayer, int *nset,
- HC_CPREC *zlabel,int *ivec,
- FILE *in, hc_boolean short_format,
- hc_boolean binary,hc_boolean verbose)
+ HC_CPREC *zlabel,
+ int *ivec,
+ FILE *in,
+ hc_boolean short_format,
+ hc_boolean binary,
+ hc_boolean verbose)
{
int input1[2],input2[3];
HC_BIN_PREC fz;
@@ -1246,7 +1250,7 @@
/* print lon lat z[i] */
fprintf(out,"%11g %11g %11g\t",lon,lat,z);
}
- for(k=0;k<shps;k++) /* loop through all scalars */
+ for(k=0;k < shps;k++) /* loop through all scalars */
fprintf(out,"%11g ",data[j+exp[0].npoints*k]);
fprintf(out,"\n");
} /* end points in lateral space loop */
Modified: mc/1D/hc/trunk/sh_rick.h
===================================================================
--- mc/1D/hc/trunk/sh_rick.h 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/sh_rick.h 2007-08-27 04:43:01 UTC (rev 7893)
@@ -68,7 +68,7 @@
int nlat,nlon,lmsize,lmsize2,nlonm1;
// logic flags
my_boolean initialized,computed_legendre,
- vector_sh_fac_init;
+ vector_sh_fac_init,sin_cos_saved;
// init
my_boolean was_called;
Modified: mc/1D/hc/trunk/sh_syn.c
===================================================================
--- mc/1D/hc/trunk/sh_syn.c 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/sh_syn.c 2007-08-27 04:43:01 UTC (rev 7893)
@@ -21,11 +21,11 @@
switches
*/
hc_boolean verbose = TRUE, short_format = FALSE ,short_format_ivec = FALSE ,binary = FALSE;
- int regular_basis = 0;
+ double regular_basis = 0;
/* */
- float *data,*theta,*phi,x,y;
+ float *data,*theta,*phi;
/* spacing for irregular output */
- static float dphi = 2.0;
+ double dphi,x,y;
HC_PREC fac[3] = {1.,1.,1.},zlabel;
SH_RICK_PREC *dummy;
struct sh_lms *exp;
@@ -44,15 +44,24 @@
short_format_ivec = TRUE;
}
if(argc > 3)
- sscanf(argv[3],"%i",®ular_basis);
+ sscanf(argv[3],"%lf",®ular_basis);
if((argc > 4)||(argc<0)){
- fprintf(stderr,"usage: %s [short_format, %i] [short_ivec,%i] [regular_basis,%i]\n",
+ fprintf(stderr,"usage: %s [short_format, %i] [short_ivec,%i] [regular_basis,%g]\n",
argv[0],short_format,short_format_ivec,regular_basis);
+ fprintf(stderr,"short_format:\n\t0: expects regular format with long header\n");
+ fprintf(stderr,"\t1: expects short format with only lmax in header\n\n");
+ fprintf(stderr,"short_ivec:\n\t0: for short format, expect AB for scalar expansion\n");
+ fprintf(stderr,"\t1: for short format, expect poloidal toroidal AP BP AT BT for vector expansion\n\n");
+ fprintf(stderr,"regular_basis:\n\t0: use Gauss latitudes and FFT divided longitudes dependening on lmax\n");
+ fprintf(stderr,"\t>0 use even spacing with regular_basis degree inrecrments on the globe\n\n");
+ fprintf(stderr,"The output format will depend on the type of SH input.\n");
+ fprintf(stderr,"\tfor scalara: lon lat scalar if a single SH is read in, else lon lat zlabel scalar.\n");
+ fprintf(stderr,"\tfor vectors: lon lat v_theta v_phi if a single SH is read in, else lon lat zlabel v_theta v_phi.\n\n\n");
exit(-1);
}
if(verbose)
- fprintf(stderr,"%s: waiting to read spherical harmonic coefficients from stdin, use -h for help\n",
+ fprintf(stderr,"%s: waiting to read spherical harmonic coefficients from stdin\n",
argv[0]);
while(sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer,&nset,
&zlabel,&ivec,stdin,short_format,
@@ -66,41 +75,37 @@
argv[0],lmax,ivec,zlabel);
/* input and init */
- sh_allocate_and_init(&exp,shps,lmax,type,ivec,verbose,regular_basis);
+ sh_allocate_and_init(&exp,shps,lmax,type,ivec,verbose,((regular_basis>0)?(1):(0)));
sh_read_coefficients_from_file(exp,shps,-1,stdin,binary,fac,verbose);
- if(regular_basis){
+ if(regular_basis > 0){
/*
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;
- }
-
+ dphi = regular_basis;
+ if(verbose)
+ fprintf(stderr,"sh_syn: using regular spaced grid with %g deg spacing\n",dphi);
+ /* */
+ dphi = DEG2RAD(dphi);
+ nphi = (TWOPI-dphi)/dphi + 1;
+ ntheta = nphi/2;
+ 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;
hc_svecalloc(&data,npoints * shps,"sh_shsyn data");
+ /* compute the expansion */
sh_compute_spatial_irreg(exp,ivec,FALSE,&dummy,
- theta,ntheta,phi,nphi,data,verbose,FALSE);
+ 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);
+ sh_print_irreg_spatial_data_to_file(exp,shps,data,
+ (nset>1)?(TRUE):(FALSE),
+ zlabel, theta,ntheta,
+ phi,nphi,stdout);
@@ -109,7 +114,9 @@
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);
+ sh_print_spatial_data_to_file(exp,shps,data,
+ (nset>1)?(TRUE):(FALSE),
+ zlabel,stdout);
}
Modified: mc/1D/hc/trunk/vdepth.dat
===================================================================
--- mc/1D/hc/trunk/vdepth.dat 2007-08-24 21:55:23 UTC (rev 7892)
+++ mc/1D/hc/trunk/vdepth.dat 2007-08-27 04:43:01 UTC (rev 7893)
@@ -16,20 +16,7 @@
789.525
645.975
502.425
-450
-410
-375
358.875
-325
-300
-275
-250
215.325
-200
-175
-150
-125
-100.5
71.775
-50
0
More information about the cig-commits
mailing list