[cig-commits] r12020 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Sat May 24 19:03:21 PDT 2008
Author: becker
Date: 2008-05-24 19:03:20 -0700 (Sat, 24 May 2008)
New Revision: 12020
Modified:
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/hc_misc.c
mc/1D/hc/trunk/hc_polsol.c
mc/1D/hc/trunk/itg-hc-geoid.ab
mc/1D/hc/trunk/main.c
mc/1D/hc/trunk/make_tar
mc/1D/hc/trunk/rick_sh_c.c
mc/1D/hc/trunk/sh_exp.c
Log:
Added basic viscosity parameter space exploration for four layer
geoid inversions.
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2008-05-24 02:57:03 UTC (rev 12019)
+++ mc/1D/hc/trunk/hc.h 2008-05-25 02:03:20 UTC (rev 12020)
@@ -151,11 +151,21 @@
hc_boolean sol_binary_out; /* binary or ASCII output of SH expansion */
hc_boolean print_spatial; /* print the spatial solution */
hc_boolean compute_geoid; /* compute and print the geoid */
+ hc_boolean compute_geoid_correlations; /* compute correlations only */
int solution_mode; /* velocity or stress */
+ hc_boolean solver_mode;
+ hc_boolean visc_init_mode;
+ HC_PREC elayer[4];
+
hc_boolean read_short_dens_sh; /* short SH format for density
files? */
- hc_boolean read_dens_scale_from_file; /* read the density/velocity scaling from file? */
+ hc_boolean dd_dens_scale; /* read the density/velocity scaling from file? */
+ HC_PREC *rdf,*sdf;
+ int ndf;
+ struct sh_lms *ref_geoid;
+ HC_PREC rlayer[3]; /* for four layer approaches (first
+ layer radius is CMB radius) */
hc_boolean print_pt_sol; /* output of p[6] and t[2] vectors */
char visc_filename[HC_CHAR_LENGTH]; /* name of viscosity profile file */
@@ -163,6 +173,7 @@
char dens_filename[HC_CHAR_LENGTH]; /* name of density model file */
char prem_model_filename[HC_CHAR_LENGTH]; /* PREM model filename */
char dens_scaling_filename[HC_CHAR_LENGTH]; /* */
+ char ref_geoid_file[HC_CHAR_LENGTH]; /* reference geoid */
};
/*
@@ -209,6 +220,7 @@
need to call with
dens_fac_changed == TRUE)
*/
+ struct sh_lms *dens_anom_orig;
HC_PREC dens_scale; /* scale for density file */
@@ -262,6 +274,7 @@
toroidal solution
*/
hc_boolean initialized,const_init,visc_init,dens_init,pvel_init; /* logic flags */
+ hc_boolean orig_danom_saved;
/* sqrt(l(l+1)) and 1/lfac factors */
HC_PREC *lfac,*ilfac;
int lfac_init;
@@ -310,9 +323,20 @@
init and assignment modes
*/
-#define HC_INIT_FROM_FILE 0
+#define HC_INIT_E_FROM_FILE 0 /* viscosity */
+#define HC_INIT_E_FOUR_LAYERS 1
+#define HC_INIT_D_FROM_FILE 0 /* density */
+#define HC_RESCALE_D 1
+#define HC_INIT_DD_FROM_FILE 0 /* depth dependent density scaling */
+#define HC_INIT_DD_FOUR_LAYERS 1 /* */
+
+#define HC_INIT_P_FROM_FILE 0 /* plates */
+
+#define HC_SOLVER_MODE_DEFAULT 0
+#define HC_SOLVER_MODE_VISC_SCAN 1
+
/*
declarations
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2008-05-24 02:57:03 UTC (rev 12019)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2008-05-25 02:03:20 UTC (rev 12020)
@@ -40,14 +40,16 @@
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 *);
-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, unsigned short, unsigned short, char *);
+void hc_assign_viscosity(struct hcs *, int, double [4], struct hc_parameters *);
+void hc_assign_density(struct hcs *, unsigned short, int, char *, int, unsigned short, unsigned short, unsigned short, unsigned short, unsigned short, int, double *, double *, unsigned short);
double hc_find_dens_scale(double, double, unsigned short, double *, double *, int);
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);
void hc_get_blank_expansions(struct sh_lms **, int, int, char *);
void hc_struc_free(struct hcs **);
+void hc_assign_dd_scaling(int, double [4], struct hc_parameters *, double);
+void hc_read_geoid(struct hc_parameters *);
/* hc_input.c */
void hc_read_sh_solution(struct hcs *, struct sh_lms **, FILE *, unsigned short, unsigned short);
/* hc_matrix.c */
@@ -73,6 +75,7 @@
char *hc_name_boolean(unsigned short);
unsigned short hc_toggle_boolean(unsigned short *);
void hc_advance_argument(int *, int, char **);
+void hc_compute_correlation(struct sh_lms *, struct sh_lms *, double *, int, unsigned short);
/* hc_output.c */
void hc_print_spectral_solution(struct hcs *, struct sh_lms *, FILE *, int, unsigned short, unsigned short);
void hc_print_sh_scalar_field(struct sh_lms *, FILE *, unsigned short, unsigned short, unsigned short);
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2008-05-24 02:57:03 UTC (rev 12019)
+++ mc/1D/hc/trunk/hc_init.c 2008-05-25 02:03:20 UTC (rev 12020)
@@ -24,22 +24,41 @@
p->free_slip = TRUE; /* free slip? */
p->no_slip = FALSE; /* no slip boundary condition? */
p->platebc = FALSE; /* plate velocities? */
-
p->compute_geoid = TRUE; /* compute the geoid? */
+ p->compute_geoid_correlations = FALSE; /* compute the geoid
+ correlation with
+ refernece only */
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->read_dens_scale_from_file = FALSE; /* read a depth-dependent scaling from file? */
+
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 */
+ p->solver_mode = HC_SOLVER_MODE_DEFAULT ;
+
p->print_pt_sol = FALSE;
p->print_spatial = FALSE; /* by default, only print the spectral solution */
+ /* for four layer approaches */
+ p->rlayer[0] = HC_ND_RADIUS(660);
+ p->rlayer[1] = HC_ND_RADIUS(410);
+ p->rlayer[2] = HC_ND_RADIUS(100);
+ /* default viscosities for the four layers, in units of reference viscosity*/
+ p->elayer[0] = 50.; p->elayer[1] = 1.; p->elayer[2] = 0.1; p->elayer[3] = 50.;
+ p->visc_init_mode = HC_INIT_E_FROM_FILE; /* by default, read viscosity from file */
+
+
/*
+ depth dependent scaling of density files?
+ */
+ p->dd_dens_scale = FALSE; /* read a depth-dependent scaling from file? */
+ p->ndf = 0;
+ p->rdf = p->sdf = NULL;
+ /*
filenames
@@ -62,7 +81,7 @@
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)->dens_init = (*hc)->pvel_init = (*hc)->orig_danom_saved = FALSE;
hc_init_polsol_struct(&((*hc)->psp));
/*
assign NULL pointers to allow reallocating
@@ -102,6 +121,7 @@
struct hc_parameters *p)
{
int dummy=0;
+ HC_PREC dd_dummy[4]={1,1,1,1};
/* mechanical boundary condition */
if(p->free_slip){
if(p->no_slip || p->platebc)
@@ -116,13 +136,23 @@
/*
start by reading in physical constants and PREM model
*/
- hc_init_constants(hc,p->dens_anom_scale,
- p->prem_model_filename,p->verbose);
+ hc_init_constants(hc,p->dens_anom_scale,p->prem_model_filename,p->verbose);
/*
initialize viscosity structure from file
*/
- hc_assign_viscosity(hc,HC_INIT_FROM_FILE,p->visc_filename,p->verbose);
+ hc_assign_viscosity(hc,p->visc_init_mode,p->elayer,p);
+ /*
+ initialize possible depth dependent scaling of density model
+ */
+ hc_assign_dd_scaling(HC_INIT_DD_FROM_FILE,dd_dummy,p,hc->r_cmb);
+
+
+ if(!p->dd_dens_scale)
+ if(p->verbose)
+ fprintf(stderr,"hc_init_main: using constant dln\\rho/dln input density scaling of %g\n",
+ hc->dens_scale);
+
if(p->no_slip && (!p->platebc)){
/*
@@ -137,13 +167,14 @@
/*
read in the densities first to determine L from the density expansion
*/
- hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,
+ hc_assign_density(hc,p->compressible,HC_INIT_D_FROM_FILE,
p->dens_filename,-1,FALSE,FALSE,p->verbose,p->read_short_dens_sh,
- p->read_dens_scale_from_file,p->dens_scaling_filename);
+ p->dd_dens_scale,p->ndf,p->rdf,p->sdf,
+ (p->solver_mode == HC_SOLVER_MODE_VISC_SCAN)?(TRUE):(FALSE));
/*
assign all zeroes up to the lmax of the density expansion
*/
- hc_assign_plate_velocities(hc,HC_INIT_FROM_FILE,p->pvel_filename,TRUE,hc->dens_anom[0].lmax,
+ hc_assign_plate_velocities(hc,HC_INIT_P_FROM_FILE,p->pvel_filename,TRUE,hc->dens_anom[0].lmax,
FALSE,p->verbose);
}else if(p->platebc){
/*
@@ -158,11 +189,11 @@
/*
read in velocities, which will determine the solution lmax
*/
- hc_assign_plate_velocities(hc,HC_INIT_FROM_FILE,p->pvel_filename,FALSE,dummy,FALSE,p->verbose);
+ hc_assign_plate_velocities(hc,HC_INIT_P_FROM_FILE,p->pvel_filename,FALSE,dummy,FALSE,p->verbose);
/* then read in the density anomalies */
- hc_assign_density(hc,p->compressible,HC_INIT_FROM_FILE,p->dens_filename,hc->pvel[0].lmax,
- FALSE,FALSE,p->verbose,p->read_short_dens_sh,
- p->read_dens_scale_from_file,p->dens_scaling_filename);
+ hc_assign_density(hc,p->compressible,HC_INIT_D_FROM_FILE,p->dens_filename,hc->pvel[0].lmax,
+ FALSE,FALSE,p->verbose,p->read_short_dens_sh, p->dd_dens_scale,p->ndf,p->rdf,p->sdf,
+ (p->solver_mode == HC_SOLVER_MODE_VISC_SCAN)?(TRUE):(FALSE));
}else if(p->free_slip){
/*
@@ -174,9 +205,9 @@
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->read_short_dens_sh,
- p->read_dens_scale_from_file,p->dens_scaling_filename);
+ hc_assign_density(hc,p->compressible,HC_INIT_D_FROM_FILE,p->dens_filename,-1,FALSE,FALSE,
+ p->verbose,p->read_short_dens_sh, p->dd_dens_scale,p->ndf,p->rdf,p->sdf,
+ (p->solver_mode == HC_SOLVER_MODE_VISC_SCAN)?(TRUE):(FALSE));
}else{
HC_ERROR("hc_init","boundary condition logic error");
}
@@ -308,14 +339,14 @@
p->dens_anom_scale);
fprintf(stderr,"-dsf\tfile\tread depth dependent density scaling from file\n");
fprintf(stderr,"\t\t(overrides -ds, %s), use pdens.py to edit\n\n",
- hc_name_boolean(p->read_dens_scale_from_file));
+ hc_name_boolean(p->dd_dens_scale));
fprintf(stderr,"Earth model options:\n");
fprintf(stderr,"-prem\tname\tset Earth model to name (%s)\n",
p->prem_model_filename);
fprintf(stderr,"-vf\tname\tviscosity structure filename (%s), use pvisc.py to edit\n",
p->visc_filename);
- fprintf(stderr,"\t\tThis file is in non_dim_radius viscosity[Pas] format\n\n");
-
+ fprintf(stderr,"\t\tThis file is in non_dim_radius viscosity[Pas] format\n");
+ fprintf(stderr,"\t\tif name is \"visc_scan\", will loop through a four layer viscosity scan\n\n");
fprintf(stderr,"boundary condition options:\n");
fprintf(stderr,"-fs\t\tperform free slip computation (%s)\n",hc_name_boolean(p->free_slip));
fprintf(stderr,"-ns\t\tperform no slip computation (%s)\n",hc_name_boolean(p->no_slip));
@@ -326,7 +357,10 @@
fprintf(stderr,"solution procedure and I/O options:\n");
fprintf(stderr,"-ng\t\tdo not compute and print the geoid (%s)\n",
hc_name_boolean(!p->compute_geoid));
- fprintf(stderr,"-pptsol\t\tprint pol[6] and tor[2] solution vectors (%s)\n",
+ fprintf(stderr,"-rg\tname\tcompute correlation of computed geoid with that in file \"name\",\n");
+ fprintf(stderr,"\t\tthis will not print out the geoid file, but only correlations (%s)\n",
+ hc_name_boolean(p->compute_geoid_correlations));
+ fprintf(stderr,"-pptsol\t\tprint pol[6] and tor[2] solution vectors (%s)\n",
hc_name_boolean(p->print_pt_sol));
fprintf(stderr,"-px\t\tprint the spatial solution to file (%s)\n",
hc_name_boolean(p->print_spatial));
@@ -354,9 +388,20 @@
hc_toggle_boolean(&p->print_pt_sol);
}else if(strcmp(argv[i],"-ng")==0){ /* do not compute geoid */
hc_toggle_boolean(&p->compute_geoid);
+ }else if(strcmp(argv[i],"-rg")==0){ /* compute geoid correlations */
+ hc_toggle_boolean(&p->compute_geoid_correlations);
+ p->compute_geoid = TRUE;
+ hc_advance_argument(&i,argc,argv); /* filename */
+ strncpy(p->ref_geoid_file,argv[i],HC_CHAR_LENGTH);
+ hc_read_geoid(p);
}else if(strcmp(argv[i],"-vf")==0){ /* viscosity filename */
hc_advance_argument(&i,argc,argv);
strncpy(p->visc_filename,argv[i],HC_CHAR_LENGTH);
+ if(strcmp(p->visc_filename,"visc_scan")==0)
+ { /* inversion mode */
+ p->solver_mode = HC_SOLVER_MODE_VISC_SCAN;
+ p->visc_init_mode = HC_INIT_E_FOUR_LAYERS;
+ }
}else if(strcmp(argv[i],"-prem")==0){ /* PREM filename */
hc_advance_argument(&i,argc,argv);
strncpy(p->prem_model_filename,argv[i],HC_CHAR_LENGTH);
@@ -368,7 +413,7 @@
hc_advance_argument(&i,argc,argv);
strncpy(p->dens_filename,argv[i],HC_CHAR_LENGTH);
}else if(strcmp(argv[i],"-dsf")==0){ /* density scaling filename */
- hc_toggle_boolean(&p->read_dens_scale_from_file);
+ hc_toggle_boolean(&p->dd_dens_scale);
hc_advance_argument(&i,argc,argv);
strncpy(p->dens_scaling_filename,argv[i],HC_CHAR_LENGTH);
}else if(strcmp(argv[i],"-visc")==0){ /* viscosity filename */
@@ -408,14 +453,30 @@
ascending
*/
-void hc_assign_viscosity(struct hcs *hc,int mode,char filename[HC_CHAR_LENGTH],
- hc_boolean verbose)
+void hc_assign_viscosity(struct hcs *hc,int mode,
+ HC_PREC elayer[4],struct hc_parameters *p)
{
FILE *in;
+ int i;
char fstring[100];
HC_PREC mean,mweight,rold,mws;
switch(mode){
- case HC_INIT_FROM_FILE:
+ case HC_INIT_E_FOUR_LAYERS:
+ /* initialize a four layer viscosity structure, viscosity values
+ for 2871-660, 660-410, 410-100,100-0 should be given in units
+ of visnor [1e21] as elayer[4]
+ */
+ hc_vecrealloc(&hc->rvisc,4,"hc_assign_viscosity");
+ hc_vecrealloc(&hc->visc,4,"hc_assign_viscosity");
+ hc->nvis = 4;
+ hc->rvisc[0] = hc->r_cmb;hc->rvisc[1] = p->rlayer[0];hc->rvisc[2] = p->rlayer[1];hc->rvisc[3] = p->rlayer[2];
+ for(i=0;i < hc->nvis;i++)
+ hc->visc[i] = elayer[i];
+ if(p->verbose)
+ fprintf(stderr,"hc_assign_viscosity: assigned four layer viscosity: %.2e %.2e %.2e %.2e\n",
+ hc->visc[0],hc->visc[1],hc->visc[2],hc->visc[3]);
+ break;
+ case HC_INIT_E_FROM_FILE:
/*
init from file part
@@ -432,7 +493,7 @@
from bottom to top
*/
- in = ggrd_open(filename,"r","hc_assign_viscosity");
+ in = ggrd_open(p->visc_filename,"r","hc_assign_viscosity");
hc_vecrealloc(&hc->rvisc,1,"hc_assign_viscosity");
hc_vecrealloc(&hc->visc,1,"hc_assign_viscosity");
hc->nvis = 0;mean = 0.0;mws = 0.0;
@@ -452,7 +513,7 @@
hc->rvisc[hc->nvis], hc->r_cmb);
exit(-1);
}
- if(verbose){
+ if(p->verbose){
/* weighted mean, should use volume, really, but this is just
for information */
mweight = ( hc->rvisc[hc->nvis] - rold);
@@ -477,7 +538,7 @@
hc->rvisc[hc->nvis-1]);
exit(-1);
}
- if(verbose){
+ if(p->verbose){
/* last entry */
mweight = ( 1.0 - hc->rvisc[hc->nvis-1]);
mws += mweight;
@@ -486,7 +547,7 @@
mean = exp(mean/mws);
fprintf(stderr,"hc_assign_viscosity: read %i layered viscosity[Pas] from %s\n",
- hc->nvis,filename);
+ hc->nvis,p->visc_filename);
fprintf(stderr,"hc_assign_viscosity: rough estimate of mean viscosity: %g x %g = %g Pas\n",
mean, hc->visnor, mean*hc->visnor);
}
@@ -524,21 +585,21 @@
*/
void hc_assign_density(struct hcs *hc,
- hc_boolean compressible,int mode,
- char *filename,int nominal_lmax,
+ hc_boolean compressible,int mode, /* compressible computation? assignment mode? */
+ char *filename,int nominal_lmax, /* input density file name */
hc_boolean layer_structure_changed,
hc_boolean density_in_binary,
hc_boolean verbose,
hc_boolean use_short_format,
- hc_boolean read_dens_scale_from_file,
- char *dens_scaling_filename)
+ hc_boolean dd_dens_scale, /* depth dependent scaling ? */
+ int ndf,HC_PREC *rdf,HC_PREC *sdf, /* depth dependent scaling factors */
+ hc_boolean save_orig_danom) /* save the original density anomalies for later rescaling */
{
FILE *in;
- int type,lmax,shps,ilayer,nset,ivec,i,j,ndf;
- HC_PREC *dtop,*dbot,zlabel,local_scale,dens_scale[0],rho0,*rdf=NULL,*sdf=NULL,smean;
- double dtmp[2];
+ int type,lmax,shps,ilayer,nset,ivec,i,j;
+ HC_PREC *dtop,*dbot,zlabel,local_scale,dens_scale[0],rho0;
hc_boolean reported = FALSE,read_on;
-
+ double dtmp[3];
hc->compressible = compressible;
hc->inho = 0;
if(hc->dens_init) /* clear old expansions, if
@@ -548,39 +609,17 @@
if(!hc->prem_init)
HC_ERROR("hc_assign_density","assign 1-D reference model (PREM) first");
switch(mode){
- case HC_INIT_FROM_FILE:
+ case HC_RESCALE_D:
+ /* resuse old density model, apply new scaling */
+ if(!hc->orig_danom_saved)
+ HC_ERROR("hc_assign_density","trying to rescale original density anomaly model, but it was not saved");
+ for(i=0;i<hc->inho;i++){
+ sh_aexp_equals_bexp_coeff((hc->dens_anom+i),(hc->dens_anom_orig+i));
+ }
+ break;
+ case HC_INIT_D_FROM_FILE:
if(hc->dens_init)
HC_ERROR("hc_assign_density","really read dens anomalies again from file?");
- if(read_dens_scale_from_file){
- /* depth depending scaling, format same as for viscosity file */
- if(verbose)
- fprintf(stderr,"hc_assign_density: reading depth dependent dln\\rho/dln density scaling from %s\n",
- dens_scaling_filename);
- ndf=0;smean = 0.0;
- in = ggrd_open(dens_scaling_filename,"r","hc_assign_density");
- while(fscanf(in,"%lf %lf",dtmp,(dtmp+1)) == 2){
- hc_vecrealloc(&rdf,(1+ndf),"hc_assign_density");
- hc_vecrealloc(&sdf,(1+ndf),"hc_assign_density");
- rdf[ndf] = dtmp[0];sdf[ndf] = dtmp[1];
- smean+=sdf[ndf];
- ndf++;
- }
- fclose(in);
- if(!ndf){
- fprintf(stderr,"hc_assign_density: error: did not read any density scaling factors from %s\n",
- dens_scaling_filename);
- exit(-1);
- }
- smean /= (HC_PREC)ndf;
- if(verbose)
- fprintf(stderr,"hc_assign_density: read scaling on %i layers, rough mean: %g\n",ndf,smean);
- /* end init */
- }else{
- if(verbose)
- fprintf(stderr,"hc_assign_density: using constant dln\\rho/dln input density scaling of %g\n",
- hc->dens_scale);
- ndf = 0;
- }
/*
read in density anomalies in spherical harmonics format for
@@ -641,7 +680,7 @@
}
reported = TRUE;
if(verbose >= 2)
- fprintf(stderr,"hc_assign_density: non_dim radius d\\rho/dinput %% factor PREM \\rho layer # depth[km]\n");
+ fprintf(stderr,"hc_assign_density: non_dim radius %% factor PREM \\rho layer # depth[km]\n");
}
/*
@@ -668,15 +707,12 @@
/*
general density
*/
- /* depth dependent factor? */
- local_scale = hc_find_dens_scale(hc->rden[hc->inho],hc->dens_scale,
- read_dens_scale_from_file,rdf,sdf,ndf);
- /* total scaling factor */
- dens_scale[0] = local_scale * ((HC_PREC)HC_DENSITY_SCALING ) * rho0;
+ /* scaling factor without depth dependence */
+ dens_scale[0] = ((HC_PREC)HC_DENSITY_SCALING ) * 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],local_scale,
+ fprintf(stderr,"hc_assign_density: r: %11g anom scales: %11g x %11g = %11g\t%5i out of %i, z: %11g\n",
+ hc->rden[hc->inho],
HC_DENSITY_SCALING,rho0,dens_scale[0],hc->inho+1,nset,zlabel);
}
if(hc->inho){
@@ -702,7 +738,7 @@
*/
sh_init_expansion((hc->dens_anom+hc->inho),
(nominal_lmax > lmax) ? (nominal_lmax):(lmax),
- hc->sh_type,1,verbose,FALSE);
+ hc->sh_type,0,verbose,FALSE);
/*
read parameters and scale (put possible depth dependence of
@@ -725,6 +761,33 @@
HC_ERROR("hc_assign_density","mode undefined");
break;
}
+ if(save_orig_danom && (!hc->dens_init)){
+ /* make a copy of the original density anomaly before applying
+ depth dependent scaling, only done once per run */
+ hc_get_blank_expansions(&hc->dens_anom_orig,hc->inho+1,hc->inho,
+ "hc_assign_density");
+ for(i=0;i<hc->inho;i++){
+ sh_init_expansion((hc->dens_anom_orig+i),hc->dens_anom[0].lmax,
+ hc->sh_type,0,FALSE,FALSE);
+ sh_aexp_equals_bexp_coeff((hc->dens_anom_orig+i),(hc->dens_anom+i));
+ }
+ hc->orig_danom_saved=TRUE;
+ }
+ /*
+
+ scale with possibly depth dependent factor
+
+ */
+ for(i=0;i < hc->inho;i++){
+ /* depth dependent factor? */
+ local_scale = hc_find_dens_scale(hc->rden[i],hc->dens_scale,dd_dens_scale,rdf,sdf,ndf);
+ sh_scale_expansion((hc->dens_anom+i),local_scale);
+ if(verbose >= 2){
+ fprintf(stderr,"hc_assign_density: r: %11g additional %s d\\rho/dinput: %11g \tlayer %5i out of %i\n",
+ hc->rden[i],(dd_dens_scale)?("depth-dependent"):("constant"),local_scale,i,hc->inho);
+ }
+ }
+
if((!hc->dens_init)||(layer_structure_changed)){
/*
@@ -794,11 +857,7 @@
free(dbot);free(dtop);
} /* end layer structure part */
- /* free density scaling factors, if used */
- if(read_dens_scale_from_file && ndf){
- free(rdf);
- free(sdf);
- }
+
hc->dens_init = TRUE;
}
@@ -876,7 +935,7 @@
*/
switch(mode){
- case HC_INIT_FROM_FILE:
+ case HC_INIT_P_FROM_FILE:
/*
read velocities in pol/tor expansion format from file in
units of HC_VELOCITY_FILE_FACTOR per year, short format
@@ -1044,3 +1103,103 @@
free(*hc);
}
+
+/*
+
+assign a depth dependent density scale
+
+*/
+void hc_assign_dd_scaling(int mode, HC_PREC dlayer[4],struct hc_parameters *p,
+ HC_PREC rcmb)
+{
+ HC_PREC smean;
+ double dtmp[2];
+ int i;
+ FILE *in;
+ if(p->dd_dens_scale){
+ /*
+ depth depending scaling
+ */
+ switch(mode){
+ case HC_INIT_DD_FROM_FILE:
+ /*
+
+ read from file, format same as for viscosity
+ */
+ if(p->verbose)
+ fprintf(stderr,"hc_assign_dd_scaling: reading depth dependent dln\\rho/dln density scaling from %s\n",
+ p->dens_scaling_filename);
+ p->ndf=0;smean = 0.0;
+ in = ggrd_open(p->dens_scaling_filename,"r","hc_assign_dd_scaling");
+ while(fscanf(in,"%lf %lf",dtmp,(dtmp+1)) == 2){
+ hc_vecrealloc(&p->rdf,(1+p->ndf),"hc_assign_dd_scaling");
+ hc_vecrealloc(&p->sdf,(1+p->ndf),"hc_assign_dd_scaling");
+ p->rdf[p->ndf] = dtmp[0];p->sdf[p->ndf] = dtmp[1];
+ smean+=p->sdf[p->ndf];
+ p->ndf++;
+ }
+ fclose(in);
+ if(!p->ndf){
+ fprintf(stderr,"hc_assign_dd_scaling: error: did not read any density scaling factors from %s\n",
+ p->dens_scaling_filename);
+ exit(-1);
+ }
+ smean /= (HC_PREC)p->ndf;
+ if(p->verbose)
+ fprintf(stderr,"hc_assign_dd_density: read scaling on %i layers, rough mean: %g\n",p->ndf,smean);
+ break;
+
+ case HC_INIT_DD_FOUR_LAYERS:
+ p->ndf = 4;
+ hc_vecrealloc(&p->rdf,(p->ndf),"hc_assign_dd_scaling");
+ hc_vecrealloc(&p->sdf,(p->ndf),"hc_assign_dd_scaling");
+ p->rdf[0] = rcmb;p->rdf[1] = p->rlayer[0];p->rdf[2] = p->rlayer[1];p->rdf[3] = p->rlayer[2];
+ for(i=0;i<4;i++)
+ p->sdf[i] = dlayer[i];
+ if(p->verbose)
+ fprintf(stderr,"hc_assign_dd_density: assigned four layer density scaling factors %g %g %g %g\n",
+ p->sdf[0], p->sdf[1], p->sdf[2], p->sdf[3]);
+ break;
+ }
+ /* end init */
+ }else{
+ p->ndf = 0;
+ hc_vecrealloc(&p->rdf,(1+p->ndf),"hc_assign_dd_scaling");
+ hc_vecrealloc(&p->sdf,(1+p->ndf),"hc_assign_dd_scaling");
+ }
+}
+
+/* read in and assign reference geoid */
+void hc_read_geoid(struct hc_parameters *p)
+{
+ int type,lmax,shps,ilayer,nset,ivec;
+ double zlabel;
+ FILE *in;
+ HC_PREC fac[3]={1,1,1};
+
+ in = fopen(p->ref_geoid_file,"r");
+ if(!in){
+ fprintf(stderr,"hc_read_geoid: error, could not open ref geoid \"%s\", expecting long scalar SH format\n",
+ p->ref_geoid_file);
+ exit(-1);
+ }
+
+
+ /* read in the file */
+ sh_read_parameters_from_file(&type,&lmax,&shps,&ilayer,&nset,
+ &zlabel,&ivec,in,FALSE,
+ FALSE,p->verbose);
+ if((ivec != 0)||(shps!=1)){
+ fprintf(stderr,"hc_read_geoid: error, expecting scalar, long format SH in %s\n",
+ p->ref_geoid_file);
+ exit(-1);
+ }
+ sh_allocate_and_init(&p->ref_geoid,shps,lmax,type,ivec,p->verbose,0);
+ sh_read_coefficients_from_file(p->ref_geoid,shps,-1,in,FALSE,fac,p->verbose);
+ fclose(in);
+ if(p->verbose)
+ fprintf(stderr,"hc_read_geoid: read reference geoid from %s, L=%i\n",
+ p->ref_geoid_file,p->ref_geoid->lmax);
+
+}
+
Modified: mc/1D/hc/trunk/hc_misc.c
===================================================================
--- mc/1D/hc/trunk/hc_misc.c 2008-05-24 02:57:03 UTC (rev 12019)
+++ mc/1D/hc/trunk/hc_misc.c 2008-05-25 02:03:20 UTC (rev 12020)
@@ -194,3 +194,26 @@
}
*i += 1;
}
+
+void hc_compute_correlation(struct sh_lms *g1,struct sh_lms *g2,
+ HC_PREC *c,int mode,hc_boolean verbose)
+{
+ int lmaxg;
+ lmaxg = MIN(g1->lmax,g1->lmax);
+ lmaxg = MIN(20,lmaxg);
+ switch(mode){
+ case 0:
+ if(verbose)
+ fprintf(stderr,"hc_compute_correlation: computing 1...%i\n",lmaxg);
+ c[0] = sh_correlation_per_degree(g1,g2,1,lmaxg);
+ case 1: /* 1...20 and 4..9 correlations */
+ if(verbose)
+ fprintf(stderr,"hc_compute_correlation: computing 1...%i and 4..9 correlations\n",lmaxg);
+ c[0] = sh_correlation_per_degree(g1,g2,1,lmaxg);
+ c[1] = sh_correlation_per_degree(g1,g2,4,9);
+ break;
+ default:
+ fprintf(stderr,"sh_compute_correlation: mode %i undefined\n",mode);
+ exit(-1);
+ }
+}
Modified: mc/1D/hc/trunk/hc_polsol.c
===================================================================
--- mc/1D/hc/trunk/hc_polsol.c 2008-05-24 02:57:03 UTC (rev 12019)
+++ mc/1D/hc/trunk/hc_polsol.c 2008-05-25 02:03:20 UTC (rev 12020)
@@ -522,7 +522,7 @@
*/
el = (HC_PREC)l;
- if((!save_prop_mats) || (!hc->psp.prop_mats_init)){
+ if((!save_prop_mats) || (!hc->psp.prop_mats_init)|| (viscosity_or_layer_changed)){
//
// get all propagators now, as they only depend on l
//
Modified: mc/1D/hc/trunk/itg-hc-geoid.ab
===================================================================
--- mc/1D/hc/trunk/itg-hc-geoid.ab 2008-05-24 02:57:03 UTC (rev 12019)
+++ mc/1D/hc/trunk/itg-hc-geoid.ab 2008-05-25 02:03:20 UTC (rev 12020)
@@ -1,529 +1,529 @@
-31 0 0.00000000e+00 1 1 0
+31 0 0 1 1 0
0 0
0 0
0 0
- -1.0023972e+02 0.0000000e+00
- 5.9891817e-03 -3.3284724e-02
- 5.5032242e+01 -3.1590032e+01
- 2.1593456e+01 0.0000000e+00
- -4.5807012e+01 -5.5993746e+00
- 2.0411921e+01 -1.3964699e+01
- -1.6272945e+01 -3.1907573e+01
- -1.0318153e+01 0.0000000e+00
- 1.2095655e+01 1.0683629e+01
- 7.9072799e+00 1.4945480e+01
- -2.2353624e+01 4.5335619e+00
- -4.2529837e+00 6.9665827e+00
- 1.5491948e+00 0.0000000e+00
- 1.4194953e+00 2.1289721e+00
- 1.4710811e+01 -7.2948141e+00
- 1.0193625e+01 4.8493715e+00
- -6.6625854e+00 1.1236422e+00
- -3.9437355e+00 1.5101141e+01
- -3.3829442e+00 0.0000000e+00
- 1.7127695e+00 -5.9811384e-01
- 1.0975139e+00 -8.4326473e+00
- -1.2914459e+00 -2.0195620e-01
- -1.9406880e+00 -1.0635309e+01
- 6.0272458e+00 1.2103227e+01
- 2.1365764e-01 -5.3553200e+00
- 2.0419429e+00 0.0000000e+00
- -6.3367927e+00 -2.1460307e+00
- 7.4539687e+00 2.0980007e+00
- -5.6503147e+00 4.8981653e+00
- -6.2038337e+00 -2.7987436e+00
- -3.7172245e-02 -4.0445784e-01
- -8.0944539e+00 3.4245521e+00
- -3.4007898e-02 -5.4384912e-01
- 1.1161641e+00 0.0000000e+00
- -5.2250508e-01 -1.3287200e+00
- 1.8051152e+00 1.4727212e+00
- 4.3708783e-01 1.9393372e+00
- -5.5127463e+00 1.5748442e+00
- 5.7981540e-01 -2.0124208e+00
- -1.4881598e+00 6.9698041e+00
- -1.5173098e+00 -1.6890268e+00
- -2.7979396e+00 2.7196367e+00
- 6.3208469e-01 0.0000000e+00
- -3.2069193e+00 -4.8279216e-01
- 4.8310736e-01 -7.1511279e-01
- 3.6233970e+00 1.6754301e+00
- -2.1128020e-01 4.4900216e-01
- 3.6802850e-01 1.2191248e+00
- 1.4164891e+00 5.0300063e+00
- 2.6617036e+00 2.1865551e+00
- 4.2443317e+00 -6.7801553e-02
- 1.0728765e+00 -2.1856121e+00
- 1.2031287e+00 0.0000000e+00
- -1.8896684e+00 2.9574283e+00
- -2.1203921e+00 -1.1567510e+00
- 1.5807937e-01 3.4773796e+00
- -1.9056687e+00 -1.7828076e+00
- 1.1119633e+00 1.1418402e+00
- -8.4791117e-01 -1.7995765e+00
- -1.8639208e-01 6.8785588e-02
- 9.1588991e-01 -2.0690544e+00
- -2.8284828e+00 8.5600446e-01
- 2.2658239e+00 -5.3827078e-01
- -1.1453289e+00 0.0000000e+00
- -3.5222233e-01 6.1190386e-01
- 4.5375884e-01 -2.2334377e+00
- 6.8982169e-01 3.3577086e+00
- -8.5614554e-01 -1.4385761e+00
- -8.4417388e-01 -1.1187634e+00
- -3.5290056e-02 7.7320665e-01
- -1.0500768e-01 2.0264469e+00
- -1.4216682e-01 5.5372484e-01
- 7.0099889e-01 -9.4905590e-01
- -1.1786302e+00 -4.1559002e-01
- -1.0430348e+00 1.5717731e+00
- 8.2199635e-01 0.0000000e+00
- 1.2088857e+00 9.7381114e-01
- 3.2185273e-01 7.0147079e-01
- -8.9384810e-01 -5.6540176e-01
- -1.5279464e+00 8.6590222e-02
- -6.9659384e-01 -1.7124444e-01
- 7.0707436e-02 8.7938868e-01
- 4.2980671e-01 -8.0599407e-01
- -5.8400082e-01 3.8207998e-01
- -9.4559205e-01 -5.6315255e-01
- -1.3986121e-01 6.9799932e-01
- -2.5638217e-01 1.4405756e-01
- -5.4679272e-02 -2.5039984e-01
- 9.4140827e-01 0.0000000e+00
- 1.1605279e+00 -8.7286557e-01
- 1.2478291e+00 -1.4143773e+00
- 4.8632424e-01 -2.2038000e+00
- -8.2372445e-02 -2.6510731e-01
- -1.3168257e+00 -1.5165763e+00
- -7.9060102e-01 -1.4153100e-01
- -6.7998200e-02 1.6515394e-01
- -2.2679966e-01 -2.2239237e-01
- -5.5881454e-01 -1.0352182e+00
- 9.2739192e-01 -8.3111487e-01
- 1.0043956e+00 1.0922263e-01
- -7.0641913e-01 1.9838634e+00
- 1.3806751e+00 -1.5374596e+00
- -5.1139004e-01 0.0000000e+00
- 4.2350515e-01 -6.5108287e-01
- -8.1032108e-01 -9.1441287e-02
- -8.2375264e-01 -4.4429870e-01
- 3.6137309e-02 -5.1126369e-01
- -6.6121361e-01 3.7876745e-01
- -4.3015925e-01 5.5420826e-02
- -8.4892282e-01 8.8742584e-02
- -7.8828192e-01 -3.4849513e-01
- -7.2082834e-01 -6.4214961e-01
- 8.7534215e-01 -2.9181685e-02
- -3.5300707e-01 8.8074607e-01
- 1.9092763e-01 -7.0208861e-01
- -7.2741311e-01 -1.0185179e+00
- -1.1700691e+00 -1.0865232e-01
- 4.9454888e-02 0.0000000e+00
- -2.1272680e-01 -2.3650262e-01
- -4.6316115e-01 -6.8358100e-01
- -1.2050653e+00 -3.9848473e-01
- -9.0627958e-01 1.5371019e-01
- -2.7618838e-01 -1.7192322e-01
- 7.4114269e-01 -8.2273894e-01
- -1.3457894e+00 -1.1446745e-01
- -7.2389028e-01 5.0021875e-01
- -3.0004634e-01 -8.5708340e-01
- 2.3155261e-01 3.3141911e-01
- 2.9432425e-02 -4.1783286e-01
- -7.3127190e-01 3.5211581e-01
- 6.3991226e-01 1.0322347e-01
- 1.1728119e-01 -5.5034952e-01
- 4.2963873e-01 1.0601846e-01
- -1.0626554e-01 0.0000000e+00
- -5.9073604e-01 -7.5219992e-01
- -5.5298298e-01 6.3238694e-01
- 7.6511252e-01 4.8143035e-01
- 9.2166187e-01 1.0825979e+00
- 2.7344682e-01 7.7670129e-02
- 3.1301226e-01 -8.0306187e-01
- 1.8187535e-01 1.9518321e-01
- -4.7836912e-01 1.2198521e-01
- 5.0568423e-01 8.9492010e-01
- -2.6635245e-01 2.6028328e-01
- -4.3116111e-01 7.2203809e-02
- 4.4134301e-01 1.5172437e-01
- -3.1074854e-01 -2.3630387e-02
- -4.3642485e-01 -8.7199011e-01
- 3.2530330e-01 7.3942321e-01
- -8.6402739e-01 6.6747523e-02
- 4.3286993e-01 0.0000000e+00
- 5.7222055e-01 7.1524402e-01
- -4.5349296e-01 1.5371763e-01
- -1.4231667e-01 -1.1467386e-01
- 1.4607870e-01 5.7156509e-01
- 3.6590775e-01 -1.8115635e-01
- -2.6467330e-01 -6.6416580e-01
- -5.6346591e-01 9.8937207e-02
- 8.7985021e-01 8.1991664e-02
- -7.8486124e-02 6.2361584e-01
- -8.5634884e-02 4.1460143e-01
- 3.6135299e-01 -2.5012262e-01
- 6.4796922e-01 4.6153512e-01
- -3.7238044e-01 -4.5425784e-01
- -3.2203733e-01 2.6111300e-01
- -1.2493265e-01 -1.1829011e-01
- -6.8580343e-01 8.1277571e-02
- 7.8276203e-01 4.4836157e-01
- 1.3758441e-01 0.0000000e+00
- -1.6246553e-01 8.8653578e-01
- 3.3219747e-01 2.4440002e-01
- 1.1380511e-01 1.3195802e-01
- 1.2323413e+00 -1.7748956e-02
- -1.3482940e-01 -5.8971213e-01
- 3.0630843e-01 -2.9879688e-01
- -1.5325979e-01 -1.6847513e-01
- 6.8805975e-01 9.8009025e-02
- 4.4148981e-01 -8.1510375e-01
- 1.1765627e-01 -9.5669383e-02
- 1.5535666e-01 -4.7810298e-02
- -6.7113371e-01 -3.7361286e-01
- 1.4108442e-01 7.8822363e-01
- -1.8713371e-01 -2.8951266e-01
- 9.1297586e-01 4.5752924e-01
- 2.2932837e-01 1.4673458e-01
- -7.8581491e-02 -9.8731282e-02
- 6.7469730e-02 -2.4497702e-01
- -7.4518371e-02 0.0000000e+00
- 2.0239779e-01 -2.6925604e-02
- 8.0626376e-01 -5.3470531e-02
- 1.7056005e-01 -2.4881329e-02
- 3.5663872e-01 -1.8370028e-01
- -2.3437283e-01 -6.1928083e-01
- -1.0827940e-01 4.2314167e-01
- -1.2726791e-01 1.9721659e-01
- 6.7429486e-01 -2.2542201e-01
- -7.4179326e-02 -1.6328539e-01
- -7.6429640e-01 -1.7098240e-01
- -3.6782967e-01 -2.3365219e-01
- -5.5041660e-02 2.1333150e-01
- 1.7013037e-01 6.4349257e-01
- -1.0701688e-01 -2.8990524e-01
- 3.9682170e-01 3.1745312e-01
- -4.9056727e-01 -1.6012297e-01
- -6.4982502e-01 3.4612830e-01
- 7.8979754e-01 -2.1859958e-01
- 6.1006751e-02 -1.1725524e-01
- 4.8637174e-01 0.0000000e+00
- -1.2560648e-01 -1.5857456e-01
- 4.5778188e-01 3.8768334e-01
- 1.0717464e-01 -8.7757777e-01
- 9.6331499e-02 -5.1043398e-01
- 2.2843510e-01 1.8634291e-01
- 2.7437423e-01 -9.8451805e-02
- 4.7886297e-01 1.5912056e-02
- 1.1494619e-01 4.8589158e-02
- -3.8773739e-01 1.5820679e-01
- -7.2858971e-01 -1.0805600e-01
- -3.2598869e-01 4.3204954e-01
- -1.4590613e-01 4.0912136e-01
- -6.1842136e-01 -1.5259148e-01
- 2.6003127e-01 -3.2445687e-01
- 5.8183465e-01 1.9514898e-02
- -2.8028906e-01 -7.7480278e-03
- -1.0159805e-01 3.0979900e-01
- 3.4644586e-01 -2.0246371e-02
- 6.8362682e-02 -2.4641032e-01
- 8.4257277e-02 -2.8639602e-01
- 1.4099411e-01 0.0000000e+00
- 3.6453370e-01 -6.4729359e-01
- -1.2671197e-01 9.4276753e-02
- -4.3193234e-01 -5.1263705e-01
- -1.1527476e-01 4.2829635e-01
- -5.8795953e-02 -2.4751975e-02
- -3.0567333e-01 3.0937361e-04
- 1.4913082e-01 -5.6878143e-02
- -3.8873542e-01 4.8476479e-02
- -3.5771804e-01 -2.0480807e-01
- -2.5717672e-01 -2.5728950e-02
- -1.5553343e-01 8.0264664e-01
- -7.5039715e-02 3.3332000e-01
- 4.3683844e-01 -3.1443968e-01
- 4.5717479e-01 1.6590296e-01
- -3.9644712e-01 -2.3392280e-01
- 1.6903425e-01 -1.5010169e-01
- 1.6845145e-01 1.6643616e-01
- 5.9717756e-01 -2.4978579e-01
- 6.1368915e-01 -3.7348322e-01
- -6.0798013e-01 3.5912864e-01
- -1.9083341e-01 8.2414149e-02
- -2.4313357e-01 0.0000000e+00
- -3.5376183e-01 8.6643561e-02
- -5.9759507e-01 -2.6274564e-02
- -2.5591680e-01 -2.2833821e-01
- -8.1554047e-02 4.2248959e-01
- -2.1840226e-03 7.1852834e-03
- 2.2931683e-01 -1.2425161e-01
- -3.9399795e-01 -1.0710242e-01
- -5.3013037e-01 8.6055983e-02
- -1.5076478e-01 -1.8897630e-01
- 1.2695265e-01 5.0424418e-01
- 1.0863254e-01 3.9623281e-01
- 5.4559361e-02 -1.8517891e-01
- 3.8710795e-01 -4.4020532e-01
- 2.4523594e-01 1.8784929e-01
- -5.8128768e-01 -1.0628943e-01
- 2.5451964e-03 -1.6443931e-01
- -1.9326983e-01 3.2937547e-01
- 2.3003930e-01 -3.6558349e-01
- -3.1712794e-01 7.5853495e-02
- -3.7839421e-01 4.4178119e-01
- 5.7112537e-01 -5.4074045e-01
- -2.2738773e-01 5.2678635e-02
- -5.0216558e-01 0.0000000e+00
- -2.0697888e-01 -3.6536242e-01
- -3.2471852e-01 -1.0345367e-01
- 5.4317852e-01 3.9470194e-01
- -5.3846851e-01 1.8707528e-01
- 1.9925269e-02 -3.7148709e-03
- -2.6908585e-01 3.6390114e-01
- 1.5537910e-01 4.3179082e-02
- 1.6605406e-01 4.1138950e-03
- -4.1018077e-02 2.8989377e-01
- 3.8039655e-01 -5.1203021e-02
- -2.0825166e-01 -3.2246507e-01
- 3.7147411e-01 -2.7091279e-01
- 2.6195827e-01 1.1742585e-01
- 1.5694840e-01 -3.3020584e-02
- -4.2741988e-01 8.1016011e-02
- 1.3328128e-01 2.4988992e-01
- 1.2325442e-01 2.9178246e-01
- 1.9274155e-01 -3.3427379e-01
- 1.2352707e-01 -2.4041523e-01
- 1.7981832e-01 -1.2304166e-01
- -3.5439809e-01 -2.6488284e-01
- -4.0549333e-01 1.0679771e-01
- -6.8929364e-02 2.7199446e-01
- -1.6452126e-03 0.0000000e+00
- 6.1598456e-02 3.4576729e-02
- 2.1656760e-02 3.4188129e-01
- 1.0848300e-01 2.1962796e-01
- 1.3472063e-01 1.1480998e-01
- 1.6057970e-01 4.8337315e-01
- 8.1713394e-02 3.0691622e-02
- 1.3679490e-01 -1.1092091e-01
- 3.4758564e-01 -8.1398747e-02
- 2.6113386e-01 4.0236800e-01
- 2.4752535e-01 4.7955371e-01
- -3.4256127e-01 -4.0165714e-01
- 2.5912672e-01 -1.3938104e-01
- 6.8177105e-02 -6.6208180e-02
- -4.5022666e-01 -4.6167424e-02
- -1.4336619e-01 3.5998133e-01
- 1.9060931e-01 6.1368956e-02
- 2.6942531e-01 1.4538231e-01
- -8.3333538e-03 -2.3248882e-01
- 1.0030175e-01 1.9571124e-01
- -1.2761228e-01 1.9496716e-01
- -1.2496063e-01 -3.1853664e-01
- 9.0145493e-02 -8.8796495e-02
- 1.3966278e-01 1.9898954e-01
- 2.8557416e-01 -8.3184789e-02
- 7.2326660e-02 0.0000000e+00
- -1.4373774e-01 2.0657954e-01
- 5.1250463e-01 2.1086580e-01
- 2.1840680e-01 3.2275387e-01
- 2.2926913e-01 6.6980423e-03
- 2.4908786e-01 9.3074486e-03
- 3.6714461e-01 3.6938152e-03
- -2.1641154e-01 1.4940363e-01
- 4.6986377e-02 4.1417322e-03
- 6.8520951e-01 -5.1782089e-01
- 2.0293082e-01 -9.9635726e-02
- -4.2629007e-02 -2.3281622e-01
- -1.7507545e-01 2.6104526e-01
- -1.7966255e-01 2.6284573e-01
- -4.4944362e-01 1.4802332e-01
- 1.0029836e-01 1.6720794e-01
- 2.9453762e-02 -2.9106668e-01
- 3.4993429e-01 8.3573540e-02
- 3.4685428e-02 -3.3577240e-01
- -1.7517280e-01 -2.2123461e-01
- -1.8013856e-01 -2.1728656e-02
- -2.4124074e-01 -1.7505974e-01
- -3.1897540e-01 8.7111465e-02
- -1.9368699e-01 2.8513090e-01
- 9.5688775e-02 -1.8943624e-01
- -2.3673224e-01 -1.1162147e-01
- 1.3293139e-01 0.0000000e+00
- 3.2332651e-03 1.4999758e-01
- -3.5781366e-02 2.5952399e-01
- -3.4076948e-01 -9.1259290e-02
- 4.3037046e-01 -4.5631744e-01
- -2.9677538e-01 -1.7647508e-01
- 2.1665360e-01 -2.3722731e-01
- 3.5808119e-02 -1.0033638e-01
- 7.8301220e-02 4.0180721e-02
- 2.9542376e-01 -1.6069255e-02
- -3.3769315e-01 -1.2771729e-01
- 1.0439430e-01 -3.7797781e-02
- -3.8379533e-01 5.2594330e-02
- 3.9594115e-03 -3.4665714e-02
- 1.7298363e-01 1.8193922e-01
- 3.0137354e-01 -1.8719412e-01
- 3.5772698e-02 -1.4437099e-01
- 2.7079052e-01 -1.7859097e-01
- -2.9186725e-01 1.1047411e-01
- 5.1081924e-02 -8.5028631e-02
- 1.4058175e-01 -2.7024700e-01
- 2.0667657e-01 -4.3804324e-02
- 2.3413166e-01 1.5938631e-01
- -2.3097382e-02 -2.4193339e-01
- 2.0371618e-01 3.3349291e-01
- -8.5311078e-02 1.3101471e-02
- 1.4846744e-02 4.6851490e-02
- 8.1522432e-02 0.0000000e+00
- -6.4552348e-02 -3.5808322e-02
- 1.4891017e-02 -1.3905154e-02
- -7.9655215e-02 -2.3648980e-01
- -3.6079116e-02 1.8179148e-01
- -3.2957984e-01 -3.3818447e-01
- 1.1208214e-01 1.4465984e-01
- 2.5113196e-01 1.5578122e-01
- -1.8132246e-01 -1.9768678e-01
- -6.9175865e-02 -2.5424446e-01
- -3.0485652e-01 2.8604576e-02
- -5.0460004e-02 2.1295484e-01
- -2.4735211e-01 4.6697266e-02
- 1.0313841e-01 1.0543109e-01
- 3.4601556e-01 2.6670962e-01
- 4.8102755e-02 -3.3188471e-02
- 5.8026149e-02 7.1513738e-02
- -9.9821077e-02 -1.1387424e-02
- -4.6736487e-02 1.7771690e-01
- 6.1562968e-03 5.4146072e-02
- -2.8682434e-02 8.2577183e-02
- -1.1106656e-01 1.5008048e-01
- -1.3132786e-01 6.5216976e-02
- 1.2004280e-01 2.6168914e-01
- 8.6750924e-03 -5.0536236e-02
- -2.7286580e-01 -1.3557951e-01
- -1.5396642e-01 -5.3839884e-02
- -1.8191150e-01 -2.4222699e-02
- -2.2046753e-01 0.0000000e+00
- 1.1906588e-01 -1.9669382e-01
- -3.3860217e-01 -1.8624606e-01
- -6.5490181e-02 -2.2974317e-01
- -3.1881387e-02 2.5458028e-01
- -2.7012880e-01 1.0653731e-01
- -1.0402942e-01 2.9162905e-01
- 5.8373230e-02 -1.2751144e-01
- -6.9242677e-02 -3.9633788e-02
- -2.0939708e-01 2.8641909e-01
- -2.3975338e-01 1.7974633e-01
- 7.9958474e-02 4.1479478e-02
- -1.9567734e-02 2.3325886e-01
- -4.7966117e-02 -1.3311639e-01
- -1.8322918e-01 -2.9017098e-01
- 2.7637078e-01 4.9314586e-02
- -1.0465953e-01 -2.9959557e-01
- -3.1123194e-01 1.1242824e-01
- 1.2871703e-01 -7.5259590e-02
- -1.4151726e-01 -5.4385851e-01
- -3.4318254e-02 1.4661188e-01
- -1.5290099e-01 -1.4925315e-01
- -5.0293977e-02 -1.6491870e-01
- -1.5541011e-01 -6.0189211e-02
- 2.4953163e-01 -3.1321351e-01
- -1.6616211e-01 3.9370227e-01
- 2.7227799e-01 8.6335820e-02
- 1.8122679e-01 -2.7241928e-02
- 1.5372504e-01 1.5075162e-01
- -1.1905902e-01 0.0000000e+00
- -9.0877034e-02 3.0522796e-01
- -5.3702906e-02 -5.5087654e-02
- -8.6031208e-02 2.6117975e-01
- -5.7278062e-01 5.4140053e-02
- 1.8009665e-01 -1.3437604e-01
- 2.6282731e-01 2.3419916e-01
- 8.0985969e-02 7.1066207e-02
- -4.4169236e-01 2.1609992e-01
- -1.7402571e-03 -6.0756514e-02
- 3.1040200e-01 1.2344204e-03
- 1.4051817e-01 -1.5492372e-01
- -1.7596199e-02 -1.0956378e-01
- 3.2669770e-02 5.9492696e-02
- -6.9442552e-02 -1.0654105e-01
- 2.1396529e-01 1.5905426e-01
- 3.4814164e-03 -3.3123673e-01
- 5.1497909e-02 1.0629352e-01
- -8.0111128e-02 -6.3664864e-02
- 1.4013413e-01 -1.5378910e-01
- -1.8323311e-01 1.3052757e-01
- 2.2470267e-01 1.1595139e-01
- 2.6417959e-01 -2.7940046e-02
- 4.3123471e-02 -5.0036191e-02
- 7.0571105e-03 -5.0518762e-02
- -1.3068963e-01 -1.9722821e-01
- 1.7744805e-01 -1.5374501e-01
- 1.7444248e-01 1.5958665e-02
- 2.1509425e-01 -1.2944934e-01
- -2.8886588e-01 1.1759543e-01
- 1.4113935e-01 0.0000000e+00
- -6.1582863e-03 -2.4524286e-02
- -2.2474494e-01 -2.6472235e-02
- -8.8403096e-02 3.0988028e-01
- -3.7579355e-02 -8.4774841e-02
- 5.9392753e-02 1.2800553e-01
- -3.5905833e-03 1.8871242e-02
- -5.3980003e-02 -8.6017562e-02
- 1.0502578e-01 7.3781576e-02
- 1.8096899e-01 2.6415468e-01
- 1.1029435e-01 -1.1181426e-01
- 2.1544295e-01 -2.3435416e-01
- 3.8573832e-01 -2.4220227e-01
- -3.1118319e-01 -6.7210073e-02
- 1.0908323e-01 1.8323484e-01
- -1.4788133e-02 2.3984980e-02
- -2.3836302e-01 9.3467198e-02
- 1.5533124e-01 1.3259762e-01
- -2.5182574e-01 -1.7976165e-01
- 3.0754650e-01 -6.2032192e-02
- -1.1781631e-01 3.0094699e-01
- 2.6308474e-01 1.2508058e-01
- -1.1828106e-01 -1.9646229e-01
- -1.2165449e-01 2.4759372e-01
- -6.1759899e-02 -5.8156754e-02
- -6.9418570e-02 3.5465357e-01
- 3.0582871e-02 2.8314483e-01
- 1.7437407e-01 -2.8319778e-01
- -1.3140580e-01 -1.8257730e-01
- -8.8850098e-02 -4.2339783e-02
- 5.8401058e-02 1.9117275e-01
- 1.5249474e-01 0.0000000e+00
- -1.6064579e-01 4.1683964e-01
- 1.4125958e-01 1.1543921e-01
- 1.4879338e-01 1.7017073e-01
- 2.5697974e-01 -1.1677588e-01
- 2.5902136e-01 -1.2310657e-01
- -1.9199188e-03 -6.6651293e-03
- 3.3998968e-02 3.0275241e-02
- -7.0679608e-02 4.8273408e-02
- -2.0426367e-02 -1.3195491e-01
- 9.8577625e-02 -1.1037145e-01
- -4.4809704e-02 -4.8964954e-01
- 1.2937246e-02 9.0143222e-02
- -2.2011509e-01 -1.0528021e-01
- -1.8172858e-01 7.2961576e-02
- -9.3799749e-02 5.3325374e-02
- -1.5519934e-01 1.3860747e-01
- 6.6421023e-02 -1.5462027e-01
- -1.1933248e-02 -5.4750100e-02
- -4.8277693e-02 -1.0759573e-01
- -5.5250452e-02 1.2934224e-01
- 2.2460069e-01 -1.6764788e-01
- -2.2312880e-01 -2.5331810e-01
- -1.8191027e-01 -9.7109301e-02
- -6.5239752e-02 -8.3825634e-02
- 3.7532804e-01 3.9764376e-02
- -2.8987244e-01 4.0198286e-02
- 2.5658324e-02 -2.4896017e-01
- 2.3039758e-01 6.4204623e-02
- 4.5649074e-02 5.2468538e-02
- -1.8620312e-02 -1.8352891e-01
- 1.9910140e-01 3.5638735e-02
+ -1.0046197e+02 0.0000000e+00
+ 6.0024606e-03 -3.3358521e-02
+ 5.5154257e+01 -3.1660072e+01
+ 2.1641332e+01 0.0000000e+00
+ -4.5908573e+01 -5.6117892e+00
+ 2.0457177e+01 -1.3995660e+01
+ -1.6309025e+01 -3.1978316e+01
+ -1.0341030e+01 0.0000000e+00
+ 1.2122473e+01 1.0707316e+01
+ 7.9248115e+00 1.4978616e+01
+ -2.2403185e+01 4.5436135e+00
+ -4.2624132e+00 6.9820287e+00
+ 1.5526296e+00 0.0000000e+00
+ 1.4226425e+00 2.1336924e+00
+ 1.4743427e+01 -7.3109878e+00
+ 1.0216225e+01 4.8601233e+00
+ -6.6773573e+00 1.1261335e+00
+ -3.9524793e+00 1.5134622e+01
+ -3.3904447e+00 0.0000000e+00
+ 1.7165669e+00 -5.9943994e-01
+ 1.0999473e+00 -8.4513437e+00
+ -1.2943092e+00 -2.0240396e-01
+ -1.9449907e+00 -1.0658889e+01
+ 6.0406091e+00 1.2130062e+01
+ 2.1413135e-01 -5.3671935e+00
+ 2.0464702e+00 0.0000000e+00
+ -6.3508423e+00 -2.1507888e+00
+ 7.4704952e+00 2.1026523e+00
+ -5.6628423e+00 4.9090253e+00
+ -6.2175885e+00 -2.8049488e+00
+ -3.7254661e-02 -4.0535459e-01
+ -8.1124005e+00 3.4321448e+00
+ -3.4083299e-02 -5.4505491e-01
+ 1.1186388e+00 0.0000000e+00
+ -5.2366355e-01 -1.3316660e+00
+ 1.8091174e+00 1.4759864e+00
+ 4.3805692e-01 1.9436370e+00
+ -5.5249689e+00 1.5783359e+00
+ 5.8110093e-01 -2.0168827e+00
+ -1.4914593e+00 6.9852572e+00
+ -1.5206739e+00 -1.6927716e+00
+ -2.8041431e+00 2.7256666e+00
+ 6.3348611e-01 0.0000000e+00
+ -3.2140295e+00 -4.8386258e-01
+ 4.8417848e-01 -7.1669830e-01
+ 3.6314306e+00 1.6791448e+00
+ -2.1174864e-01 4.4999766e-01
+ 3.6884447e-01 1.2218277e+00
+ 1.4196297e+00 5.0411585e+00
+ 2.6676050e+00 2.1914030e+00
+ 4.2537420e+00 -6.7951880e-02
+ 1.0752552e+00 -2.1904579e+00
+ 1.2057962e+00 0.0000000e+00
+ -1.8938581e+00 2.9639854e+00
+ -2.1250933e+00 -1.1593157e+00
+ 1.5842986e-01 3.4850895e+00
+ -1.9098938e+00 -1.7867603e+00
+ 1.1144287e+00 1.1443718e+00
+ -8.4979112e-01 -1.8035664e+00
+ -1.8680534e-01 6.8938096e-02
+ 9.1792058e-01 -2.0736418e+00
+ -2.8347539e+00 8.5790235e-01
+ 2.2708476e+00 -5.3946421e-01
+ -1.1478683e+00 0.0000000e+00
+ -3.5300326e-01 6.1326054e-01
+ 4.5476489e-01 -2.2383896e+00
+ 6.9135113e-01 3.3651531e+00
+ -8.5804374e-01 -1.4417656e+00
+ -8.4604554e-01 -1.1212439e+00
+ -3.5368299e-02 7.7492096e-01
+ -1.0524050e-01 2.0309399e+00
+ -1.4248203e-01 5.5495253e-01
+ 7.0255311e-01 -9.5116010e-01
+ -1.1812434e+00 -4.1651145e-01
+ -1.0453473e+00 1.5752579e+00
+ 8.2381883e-01 0.0000000e+00
+ 1.2115660e+00 9.7597022e-01
+ 3.2256633e-01 7.0302605e-01
+ -8.9582990e-01 -5.6665534e-01
+ -1.5313341e+00 8.6782206e-02
+ -6.9813829e-01 -1.7162411e-01
+ 7.0864204e-02 8.8133841e-01
+ 4.3075966e-01 -8.0778107e-01
+ -5.8529564e-01 3.8292711e-01
+ -9.4768857e-01 -5.6440114e-01
+ -1.4017131e-01 6.9954689e-01
+ -2.5695060e-01 1.4437696e-01
+ -5.4800504e-02 -2.5095501e-01
+ 9.4349552e-01 0.0000000e+00
+ 1.1631010e+00 -8.7480084e-01
+ 1.2505957e+00 -1.4175132e+00
+ 4.8740249e-01 -2.2086861e+00
+ -8.2555077e-02 -2.6569509e-01
+ -1.3197453e+00 -1.5199387e+00
+ -7.9235390e-01 -1.4184480e-01
+ -6.8148962e-02 1.6552011e-01
+ -2.2730250e-01 -2.2288545e-01
+ -5.6005351e-01 -1.0375134e+00
+ 9.2944808e-01 -8.3295758e-01
+ 1.0066225e+00 1.0946479e-01
+ -7.0798537e-01 1.9882619e+00
+ 1.3837363e+00 -1.5408683e+00
+ -5.1252387e-01 0.0000000e+00
+ 4.2444413e-01 -6.5252642e-01
+ -8.1211768e-01 -9.1644026e-02
+ -8.2557902e-01 -4.4528378e-01
+ 3.6217430e-02 -5.1239724e-01
+ -6.6267961e-01 3.7960724e-01
+ -4.3111297e-01 5.5543702e-02
+ -8.5080500e-01 8.8939340e-02
+ -7.9002966e-01 -3.4926779e-01
+ -7.2242653e-01 -6.4357335e-01
+ 8.7728291e-01 -2.9246385e-02
+ -3.5378974e-01 8.8269882e-01
+ 1.9135094e-01 -7.0364525e-01
+ -7.2902589e-01 -1.0207761e+00
+ -1.1726633e+00 -1.0889321e-01
+ 4.9564537e-02 0.0000000e+00
+ -2.1319845e-01 -2.3702698e-01
+ -4.6418805e-01 -6.8509660e-01
+ -1.2077371e+00 -3.9936823e-01
+ -9.0828893e-01 1.5405099e-01
+ -2.7680073e-01 -1.7230440e-01
+ 7.4278591e-01 -8.2456307e-01
+ -1.3487732e+00 -1.1472124e-01
+ -7.2549525e-01 5.0132781e-01
+ -3.0071159e-01 -8.5898369e-01
+ 2.3206599e-01 3.3215392e-01
+ 2.9497681e-02 -4.1875925e-01
+ -7.3289324e-01 3.5289650e-01
+ 6.4133104e-01 1.0345233e-01
+ 1.1754122e-01 -5.5156973e-01
+ 4.3059130e-01 1.0625352e-01
+ -1.0650115e-01 0.0000000e+00
+ -5.9204579e-01 -7.5386766e-01
+ -5.5420903e-01 6.3378904e-01
+ 7.6680889e-01 4.8249775e-01
+ 9.2370533e-01 1.0849982e+00
+ 2.7405309e-01 7.7842335e-02
+ 3.1370625e-01 -8.0484238e-01
+ 1.8227859e-01 1.9561596e-01
+ -4.7942974e-01 1.2225567e-01
+ 5.0680541e-01 8.9690427e-01
+ -2.6694299e-01 2.6086036e-01
+ -4.3211705e-01 7.2363896e-02
+ 4.4232154e-01 1.5206076e-01
+ -3.1143752e-01 -2.3682779e-02
+ -4.3739247e-01 -8.7392344e-01
+ 3.2602454e-01 7.4106262e-01
+ -8.6594307e-01 6.6895512e-02
+ 4.3382967e-01 0.0000000e+00
+ 5.7348925e-01 7.1682983e-01
+ -4.5449842e-01 1.5405844e-01
+ -1.4263221e-01 -1.1492811e-01
+ 1.4640258e-01 5.7283233e-01
+ 3.6671902e-01 -1.8155800e-01
+ -2.6526012e-01 -6.6563835e-01
+ -5.6471519e-01 9.9156565e-02
+ 8.8180097e-01 8.2173451e-02
+ -7.8660140e-02 6.2499849e-01
+ -8.5824750e-02 4.1552066e-01
+ 3.6215416e-01 -2.5067718e-01
+ 6.4940586e-01 4.6255842e-01
+ -3.7320606e-01 -4.5526500e-01
+ -3.2275133e-01 2.6169192e-01
+ -1.2520965e-01 -1.1855238e-01
+ -6.8732396e-01 8.1457776e-02
+ 7.8449753e-01 4.4935566e-01
+ 1.3788946e-01 0.0000000e+00
+ -1.6282574e-01 8.8850136e-01
+ 3.3293400e-01 2.4494189e-01
+ 1.1405743e-01 1.3225059e-01
+ 1.2350736e+00 -1.7788308e-02
+ -1.3512834e-01 -5.9101961e-01
+ 3.0698756e-01 -2.9945936e-01
+ -1.5359958e-01 -1.6884867e-01
+ 6.8958528e-01 9.8226325e-02
+ 4.4246866e-01 -8.1691095e-01
+ 1.1791713e-01 -9.5881496e-02
+ 1.5570110e-01 -4.7916300e-02
+ -6.7262171e-01 -3.7444121e-01
+ 1.4139723e-01 7.8997124e-01
+ -1.8754861e-01 -2.9015456e-01
+ 9.1500006e-01 4.5854365e-01
+ 2.2983682e-01 1.4705991e-01
+ -7.8755718e-02 -9.8950184e-02
+ 6.7619320e-02 -2.4552017e-01
+ -7.4683590e-02 0.0000000e+00
+ 2.0284654e-01 -2.6985302e-02
+ 8.0805137e-01 -5.3589083e-02
+ 1.7093820e-01 -2.4936495e-02
+ 3.5742944e-01 -1.8410757e-01
+ -2.3489247e-01 -6.2065387e-01
+ -1.0851947e-01 4.2407984e-01
+ -1.2755008e-01 1.9765385e-01
+ 6.7578987e-01 -2.2592181e-01
+ -7.4343793e-02 -1.6364742e-01
+ -7.6599096e-01 -1.7136150e-01
+ -3.6864520e-01 -2.3417023e-01
+ -5.5163696e-02 2.1380449e-01
+ 1.7050757e-01 6.4491929e-01
+ -1.0725415e-01 -2.9054800e-01
+ 3.9770151e-01 3.1815696e-01
+ -4.9165494e-01 -1.6047799e-01
+ -6.5126577e-01 3.4689572e-01
+ 7.9154864e-01 -2.1908425e-01
+ 6.1142012e-02 -1.1751521e-01
+ 4.8745010e-01 0.0000000e+00
+ -1.2588497e-01 -1.5892615e-01
+ 4.5879685e-01 3.8854289e-01
+ 1.0741226e-01 -8.7952349e-01
+ 9.6545080e-02 -5.1156569e-01
+ 2.2894158e-01 1.8675606e-01
+ 2.7498255e-01 -9.8670087e-02
+ 4.7992468e-01 1.5947336e-02
+ 1.1520104e-01 4.8696888e-02
+ -3.8859706e-01 1.5855756e-01
+ -7.3020511e-01 -1.0829557e-01
+ -3.2671146e-01 4.3300745e-01
+ -1.4622963e-01 4.1002844e-01
+ -6.1979249e-01 -1.5292980e-01
+ 2.6060780e-01 -3.2517624e-01
+ 5.8312466e-01 1.9558166e-02
+ -2.8091050e-01 -7.7652063e-03
+ -1.0182331e-01 3.1048587e-01
+ 3.4721398e-01 -2.0291260e-02
+ 6.8514252e-02 -2.4695664e-01
+ 8.4444088e-02 -2.8703100e-01
+ 1.4130671e-01 0.0000000e+00
+ 3.6534193e-01 -6.4872874e-01
+ -1.2699291e-01 9.4485779e-02
+ -4.3288999e-01 -5.1377364e-01
+ -1.1553034e-01 4.2924595e-01
+ -5.8926312e-02 -2.4806854e-02
+ -3.0635105e-01 3.1005953e-04
+ 1.4946147e-01 -5.7004250e-02
+ -3.8959731e-01 4.8583959e-02
+ -3.5851115e-01 -2.0526216e-01
+ -2.5774692e-01 -2.5785994e-02
+ -1.5587827e-01 8.0442622e-01
+ -7.5206089e-02 3.3405902e-01
+ 4.3780698e-01 -3.1513684e-01
+ 4.5818841e-01 1.6627079e-01
+ -3.9732611e-01 -2.3444144e-01
+ 1.6940902e-01 -1.5043448e-01
+ 1.6882493e-01 1.6680517e-01
+ 5.9850159e-01 -2.5033960e-01
+ 6.1504979e-01 -3.7431128e-01
+ -6.0932811e-01 3.5992488e-01
+ -1.9125652e-01 8.2596874e-02
+ -2.4367263e-01 0.0000000e+00
+ -3.5454617e-01 8.6835663e-02
+ -5.9892003e-01 -2.6332819e-02
+ -2.5648420e-01 -2.2884447e-01
+ -8.1734864e-02 4.2342631e-01
+ -2.1888649e-03 7.2012142e-03
+ 2.2982526e-01 -1.2452710e-01
+ -3.9487150e-01 -1.0733988e-01
+ -5.3130574e-01 8.6246782e-02
+ -1.5109905e-01 -1.8939529e-01
+ 1.2723412e-01 5.0536216e-01
+ 1.0887340e-01 3.9711132e-01
+ 5.4680327e-02 -1.8558948e-01
+ 3.8796623e-01 -4.4118132e-01
+ 2.4577967e-01 1.8826578e-01
+ -5.8257648e-01 -1.0652508e-01
+ 2.5508395e-03 -1.6480390e-01
+ -1.9369833e-01 3.3010575e-01
+ 2.3054934e-01 -3.6639404e-01
+ -3.1783106e-01 7.6021674e-02
+ -3.7923317e-01 4.4276069e-01
+ 5.7239164e-01 -5.4193935e-01
+ -2.2789189e-01 5.2795431e-02
+ -5.0327896e-01 0.0000000e+00
+ -2.0743778e-01 -3.6617248e-01
+ -3.2543847e-01 -1.0368304e-01
+ 5.4438282e-01 3.9557705e-01
+ -5.3966238e-01 1.8749005e-01
+ 1.9969446e-02 -3.7231073e-03
+ -2.6968245e-01 3.6470797e-01
+ 1.5572359e-01 4.3274816e-02
+ 1.6642223e-01 4.1230161e-03
+ -4.1109021e-02 2.9053651e-01
+ 3.8123994e-01 -5.1316546e-02
+ -2.0871338e-01 -3.2318002e-01
+ 3.7229772e-01 -2.7151345e-01
+ 2.6253907e-01 1.1768620e-01
+ 1.5729637e-01 -3.3093795e-02
+ -4.2836754e-01 8.1195635e-02
+ 1.3357679e-01 2.5044397e-01
+ 1.2352769e-01 2.9242938e-01
+ 1.9316889e-01 -3.3501492e-01
+ 1.2380095e-01 -2.4094826e-01
+ 1.8021701e-01 -1.2331446e-01
+ -3.5518384e-01 -2.6547013e-01
+ -4.0639237e-01 1.0703449e-01
+ -6.9082191e-02 2.7259751e-01
+ -1.6488603e-03 0.0000000e+00
+ 6.1735029e-02 3.4653391e-02
+ 2.1704776e-02 3.4263929e-01
+ 1.0872352e-01 2.2011490e-01
+ 1.3501932e-01 1.1506453e-01
+ 1.6093573e-01 4.8444486e-01
+ 8.1894565e-02 3.0759670e-02
+ 1.3709820e-01 -1.1116684e-01
+ 3.4835629e-01 -8.1579220e-02
+ 2.6171283e-01 4.0326011e-01
+ 2.4807415e-01 4.8061696e-01
+ -3.4332078e-01 -4.0254767e-01
+ 2.5970124e-01 -1.3969006e-01
+ 6.8328263e-02 -6.6354973e-02
+ -4.5122488e-01 -4.6269784e-02
+ -1.4368405e-01 3.6077946e-01
+ 1.9103192e-01 6.1505021e-02
+ 2.7002266e-01 1.4570465e-01
+ -8.3518301e-03 -2.3300428e-01
+ 1.0052413e-01 1.9614516e-01
+ -1.2789522e-01 1.9539943e-01
+ -1.2523768e-01 -3.1924289e-01
+ 9.0345359e-02 -8.8993370e-02
+ 1.3997243e-01 1.9943073e-01
+ 2.8620732e-01 -8.3369222e-02
+ 7.2487019e-02 0.0000000e+00
+ -1.4405643e-01 2.0703756e-01
+ 5.1364093e-01 2.1133332e-01
+ 2.1889104e-01 3.2346946e-01
+ 2.2977746e-01 6.7128928e-03
+ 2.4964012e-01 9.3280846e-03
+ 3.6795863e-01 3.7020049e-03
+ -2.1689136e-01 1.4973488e-01
+ 4.7090553e-02 4.1509150e-03
+ 6.8672872e-01 -5.1896898e-01
+ 2.0338075e-01 -9.9856633e-02
+ -4.2723522e-02 -2.3333240e-01
+ -1.7546362e-01 2.6162404e-01
+ -1.8006089e-01 2.6342849e-01
+ -4.5044011e-01 1.4835151e-01
+ 1.0052074e-01 1.6757867e-01
+ 2.9519066e-02 -2.9171202e-01
+ 3.5071014e-01 8.3758835e-02
+ 3.4762331e-02 -3.3651686e-01
+ -1.7556118e-01 -2.2172512e-01
+ -1.8053796e-01 -2.1776832e-02
+ -2.4177561e-01 -1.7544787e-01
+ -3.1968262e-01 8.7304604e-02
+ -1.9411642e-01 2.8576308e-01
+ 9.5900932e-02 -1.8985625e-01
+ -2.3725711e-01 -1.1186895e-01
+ 1.3322612e-01 0.0000000e+00
+ 3.2404337e-03 1.5033015e-01
+ -3.5860698e-02 2.6009940e-01
+ -3.4152502e-01 -9.1461626e-02
+ 4.3132465e-01 -4.5732916e-01
+ -2.9743337e-01 -1.7686635e-01
+ 2.1713395e-01 -2.3775328e-01
+ 3.5887511e-02 -1.0055884e-01
+ 7.8474825e-02 4.0269808e-02
+ 2.9607875e-01 -1.6104883e-02
+ -3.3844187e-01 -1.2800046e-01
+ 1.0462576e-01 -3.7881584e-02
+ -3.8464626e-01 5.2710939e-02
+ 3.9681901e-03 -3.4742573e-02
+ 1.7336716e-01 1.8234260e-01
+ 3.0204173e-01 -1.8760916e-01
+ 3.5852012e-02 -1.4469108e-01
+ 2.7139090e-01 -1.7898693e-01
+ -2.9251436e-01 1.1071905e-01
+ 5.1195180e-02 -8.5217152e-02
+ 1.4089344e-01 -2.7084618e-01
+ 2.0713480e-01 -4.3901445e-02
+ 2.3465076e-01 1.5973969e-01
+ -2.3148592e-02 -2.4246980e-01
+ 2.0416785e-01 3.3423231e-01
+ -8.5500225e-02 1.3130519e-02
+ 1.4879662e-02 4.6955367e-02
+ 8.1703180e-02 0.0000000e+00
+ -6.4695470e-02 -3.5887714e-02
+ 1.4924033e-02 -1.3935984e-02
+ -7.9831823e-02 -2.3701413e-01
+ -3.6159109e-02 1.8219454e-01
+ -3.3031056e-01 -3.3893427e-01
+ 1.1233064e-01 1.4498057e-01
+ 2.5168875e-01 1.5612661e-01
+ -1.8172447e-01 -1.9812508e-01
+ -6.9329238e-02 -2.5480816e-01
+ -3.0553243e-01 2.8667997e-02
+ -5.0571882e-02 2.1342699e-01
+ -2.4790052e-01 4.6800800e-02
+ 1.0336709e-01 1.0566485e-01
+ 3.4678273e-01 2.6730095e-01
+ 4.8209406e-02 -3.3262054e-02
+ 5.8154802e-02 7.1672295e-02
+ -1.0004239e-01 -1.1412671e-02
+ -4.6840109e-02 1.7811092e-01
+ 6.1699462e-03 5.4266122e-02
+ -2.8746028e-02 8.2760268e-02
+ -1.1131281e-01 1.5041323e-01
+ -1.3161903e-01 6.5361572e-02
+ 1.2030895e-01 2.6226935e-01
+ 8.6943264e-03 -5.0648282e-02
+ -2.7347079e-01 -1.3588011e-01
+ -1.5430779e-01 -5.3959255e-02
+ -1.8231483e-01 -2.4276404e-02
+ -2.2095634e-01 0.0000000e+00
+ 1.1932987e-01 -1.9712991e-01
+ -3.3935290e-01 -1.8665900e-01
+ -6.5635383e-02 -2.3025254e-01
+ -3.1952073e-02 2.5514473e-01
+ -2.7072772e-01 1.0677352e-01
+ -1.0426007e-01 2.9227563e-01
+ 5.8502652e-02 -1.2779416e-01
+ -6.9396199e-02 -3.9721662e-02
+ -2.0986135e-01 2.8705412e-01
+ -2.4028495e-01 1.8014485e-01
+ 8.0135754e-02 4.1571444e-02
+ -1.9611118e-02 2.3377602e-01
+ -4.8072465e-02 -1.3341153e-01
+ -1.8363542e-01 -2.9081433e-01
+ 2.7698353e-01 4.9423924e-02
+ -1.0489157e-01 -3.0025982e-01
+ -3.1192199e-01 1.1267751e-01
+ 1.2900241e-01 -7.5426452e-02
+ -1.4183103e-01 -5.4506433e-01
+ -3.4394343e-02 1.4693694e-01
+ -1.5323999e-01 -1.4958407e-01
+ -5.0405486e-02 -1.6528435e-01
+ -1.5575468e-01 -6.0322659e-02
+ 2.5008488e-01 -3.1390796e-01
+ -1.6653051e-01 3.9457516e-01
+ 2.7288167e-01 8.6527239e-02
+ 1.8162860e-01 -2.7302327e-02
+ 1.5406588e-01 1.5108586e-01
+ -1.1932299e-01 0.0000000e+00
+ -9.1078522e-02 3.0590469e-01
+ -5.3821973e-02 -5.5209791e-02
+ -8.6221952e-02 2.6175883e-01
+ -5.7405056e-01 5.4260089e-02
+ 1.8049595e-01 -1.3467397e-01
+ 2.6341004e-01 2.3471841e-01
+ 8.1165527e-02 7.1223772e-02
+ -4.4267166e-01 2.1657905e-01
+ -1.7441155e-03 -6.0891220e-02
+ 3.1109021e-01 1.2371573e-03
+ 1.4082972e-01 -1.5526721e-01
+ -1.7635212e-02 -1.0980670e-01
+ 3.2742204e-02 5.9624600e-02
+ -6.9596516e-02 -1.0677727e-01
+ 2.1443968e-01 1.5940691e-01
+ 3.4891352e-03 -3.3197113e-01
+ 5.1612088e-02 1.0652919e-01
+ -8.0288746e-02 -6.3806019e-02
+ 1.4044482e-01 -1.5413007e-01
+ -1.8363936e-01 1.3081697e-01
+ 2.2520087e-01 1.1620847e-01
+ 2.6476532e-01 -2.8001993e-02
+ 4.3219082e-02 -5.0147129e-02
+ 7.0727572e-03 -5.0630769e-02
+ -1.3097939e-01 -1.9766550e-01
+ 1.7784148e-01 -1.5408588e-01
+ 1.7482924e-01 1.5994048e-02
+ 2.1557115e-01 -1.2973634e-01
+ -2.8950633e-01 1.1785615e-01
+ 1.4145228e-01 0.0000000e+00
+ -6.1719402e-03 -2.4578660e-02
+ -2.2524324e-01 -2.6530928e-02
+ -8.8599099e-02 3.1056733e-01
+ -3.7662674e-02 -8.4962799e-02
+ 5.9524435e-02 1.2828934e-01
+ -3.5985441e-03 1.8913083e-02
+ -5.4099685e-02 -8.6208276e-02
+ 1.0525864e-01 7.3945161e-02
+ 1.8137023e-01 2.6474035e-01
+ 1.1053889e-01 -1.1206217e-01
+ 2.1592062e-01 -2.3487375e-01
+ 3.8659356e-01 -2.4273926e-01
+ -3.1187313e-01 -6.7359088e-02
+ 1.0932508e-01 1.8364109e-01
+ -1.4820920e-02 2.4038158e-02
+ -2.3889150e-01 9.3674429e-02
+ 1.5567563e-01 1.3289161e-01
+ -2.5238407e-01 -1.8016021e-01
+ 3.0822838e-01 -6.2169727e-02
+ -1.1807753e-01 3.0161424e-01
+ 2.6366804e-01 1.2535790e-01
+ -1.1854330e-01 -1.9689787e-01
+ -1.2192421e-01 2.4814268e-01
+ -6.1896829e-02 -5.8285696e-02
+ -6.9572481e-02 3.5543989e-01
+ 3.0650678e-02 2.8377261e-01
+ 1.7476068e-01 -2.8382568e-01
+ -1.3169714e-01 -1.8298210e-01
+ -8.9047092e-02 -4.2433656e-02
+ 5.8530542e-02 1.9159661e-01
+ 1.5283284e-01 0.0000000e+00
+ -1.6100197e-01 4.1776383e-01
+ 1.4157277e-01 1.1569515e-01
+ 1.4912327e-01 1.7054802e-01
+ 2.5754950e-01 -1.1703479e-01
+ 2.5959565e-01 -1.2337952e-01
+ -1.9241755e-03 -6.6799069e-03
+ 3.4074349e-02 3.0342366e-02
+ -7.0836315e-02 4.8380438e-02
+ -2.0471655e-02 -1.3224747e-01
+ 9.8796186e-02 -1.1061616e-01
+ -4.4909054e-02 -4.9073517e-01
+ 1.2965930e-02 9.0343083e-02
+ -2.2060312e-01 -1.0551363e-01
+ -1.8213150e-01 7.3123342e-02
+ -9.4007717e-02 5.3443604e-02
+ -1.5554344e-01 1.3891478e-01
+ 6.6568288e-02 -1.5496309e-01
+ -1.1959706e-02 -5.4871489e-02
+ -4.8384732e-02 -1.0783429e-01
+ -5.5372950e-02 1.2962901e-01
+ 2.2509866e-01 -1.6801958e-01
+ -2.2362351e-01 -2.5387974e-01
+ -1.8231359e-01 -9.7324607e-02
+ -6.5384398e-02 -8.4011488e-02
+ 3.7616019e-01 3.9852540e-02
+ -2.9051514e-01 4.0287411e-02
+ 2.5715213e-02 -2.4951215e-01
+ 2.3090840e-01 6.4346974e-02
+ 4.5750285e-02 5.2584868e-02
+ -1.8661596e-02 -1.8393582e-01
+ 1.9954284e-01 3.5717752e-02
Modified: mc/1D/hc/trunk/main.c
===================================================================
--- mc/1D/hc/trunk/main.c 2008-05-24 02:57:03 UTC (rev 12019)
+++ mc/1D/hc/trunk/main.c 2008-05-25 02:03:20 UTC (rev 12020)
@@ -52,12 +52,14 @@
struct hcs *model; /* main structure, make sure to initialize with
zeroes */
struct sh_lms *sol_spectral=NULL, *geoid = NULL; /* solution expansions */
- int nsol,lmax;
+ int nsol,lmax,solved;
FILE *out;
struct hc_parameters p[1]; /* parameters */
char filename[HC_CHAR_LENGTH],file_prefix[10];
float *sol_spatial = NULL; /* spatial solution,
e.g. velocities */
+ HC_PREC corr[2]; /* correlations */
+ HC_PREC vl[3][3],v1,v2,v3; /* for viscosity scans */
/*
@@ -91,7 +93,7 @@
fprintf(stderr,"%s: starting version compiled on %s\n",argv[0],__TIMESTAMP__);
#else
if(p->verbose)
- fprintf(stderr,"%s: starting\n",argv[0]);
+ fprintf(stderr,"%s: starting main program\n",argv[0]);
#endif
/*
@@ -135,75 +137,121 @@
/* make room for geoid solution */
sh_allocate_and_init(&geoid,1,model->dens_anom[0].lmax,
model->sh_type,0,p->verbose,FALSE);
-
- /*
- solve poloidal and toroidal part and sum
- */
- hc_solve(model,p->free_slip,p->solution_mode,sol_spectral,
- TRUE,TRUE,TRUE,p->print_pt_sol,p->compute_geoid,geoid,
- p->verbose);
- /*
-
- OUTPUT PART
-
- */
- /*
- output of spherical harmonics solution
- */
- switch(p->solution_mode){
- case HC_VEL:
- sprintf(file_prefix,"vel");break;
- case HC_RTRACTIONS:
- sprintf(file_prefix,"rtrac");break;
- case HC_HTRACTIONS:
- sprintf(file_prefix,"htrac");break;
- default:
- HC_ERROR(argv[0],"solution mode undefined");break;
- }
- if(p->sol_binary_out)
- sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_BINARY);
- else
- sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_ASCII);
- if(p->verbose)
- fprintf(stderr,"%s: writing spherical harmonics solution to %s\n",
- argv[0],filename);
- out = ggrd_open(filename,"w","main");
- hc_print_spectral_solution(model,sol_spectral,out,
- p->solution_mode,
- p->sol_binary_out,p->verbose);
- fclose(out);
- if(p->compute_geoid){
+ switch(p->solver_mode){
+ case HC_SOLVER_MODE_DEFAULT:
/*
- print geoid solution
+ solve poloidal and toroidal part and sum
*/
- sprintf(filename,"%s",HC_GEOID_FILE);
- if(p->verbose)
- fprintf(stderr,"%s: writing geoid to %s\n",argv[0],filename);
- out = ggrd_open(filename,"w","main");
- hc_print_sh_scalar_field(geoid,out,FALSE,FALSE,p->verbose);
- fclose(out);
- }
-
- if(p->print_spatial){
+ hc_solve(model,p->free_slip,p->solution_mode,sol_spectral,
+ TRUE,TRUE,TRUE,p->print_pt_sol,p->compute_geoid,geoid,
+ p->verbose);
/*
- we wish to use the spatial solution
-
- expand velocities to spatial base, compute spatial
- representation
-
+
+ OUTPUT PART
+
*/
- hc_compute_sol_spatial(model,sol_spectral,&sol_spatial,
- p->verbose);
/*
- output of spatial solution
+ output of spherical harmonics solution
*/
- sprintf(filename,"%s.%s",file_prefix,HC_SPATIAL_SOLOUT_FILE);
- /* print lon lat z v_r v_theta v_phi */
- hc_print_spatial_solution(model,sol_spectral,sol_spatial,
- filename,HC_LAYER_OUT_FILE,
- p->solution_mode,p->sol_binary_out,
- p->verbose);
+ switch(p->solution_mode){
+ case HC_VEL:
+ sprintf(file_prefix,"vel");break;
+ case HC_RTRACTIONS:
+ sprintf(file_prefix,"rtrac");break;
+ case HC_HTRACTIONS:
+ sprintf(file_prefix,"htrac");break;
+ default:
+ HC_ERROR(argv[0],"solution mode undefined");break;
+ }
+ if(p->sol_binary_out)
+ sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_BINARY);
+ else
+ sprintf(filename,"%s.%s",file_prefix,HC_SOLOUT_FILE_ASCII);
+ if(p->verbose)
+ fprintf(stderr,"%s: writing spherical harmonics solution to %s\n",
+ argv[0],filename);
+ out = ggrd_open(filename,"w","main");
+ hc_print_spectral_solution(model,sol_spectral,out,
+ p->solution_mode,
+ p->sol_binary_out,p->verbose);
+ fclose(out);
+ if(p->compute_geoid){
+ if(p->compute_geoid_correlations){
+ if(p->verbose)
+ fprintf(stderr,"%s: correlation for geoid with %s\n",argv[0],
+ p->ref_geoid_file);
+ hc_compute_correlation(geoid,p->ref_geoid,corr,1,p->verbose);
+ fprintf(stdout,"%10.7f %10.7f\n",corr[0],corr[1]);
+ }else{
+ /*
+ print geoid solution
+ */
+ sprintf(filename,"%s",HC_GEOID_FILE);
+ if(p->verbose)
+ fprintf(stderr,"%s: writing geoid to %s\n",argv[0],filename);
+ out = ggrd_open(filename,"w","main");
+ hc_print_sh_scalar_field(geoid,out,FALSE,FALSE,p->verbose);
+ fclose(out);
+ }
+ }
+ if(p->print_spatial){
+ /*
+ we wish to use the spatial solution
+
+ expand velocities to spatial base, compute spatial
+ representation
+
+ */
+ hc_compute_sol_spatial(model,sol_spectral,&sol_spatial,
+ p->verbose);
+ /*
+ output of spatial solution
+ */
+ sprintf(filename,"%s.%s",file_prefix,HC_SPATIAL_SOLOUT_FILE);
+ /* print lon lat z v_r v_theta v_phi */
+ hc_print_spatial_solution(model,sol_spectral,sol_spatial,
+ filename,HC_LAYER_OUT_FILE,
+ p->solution_mode,p->sol_binary_out,
+ p->verbose);
+ }
+ break; /* end default mode */
+ case HC_SOLVER_MODE_VISC_SCAN: /* scan through viscosities for a
+ four layer problem */
+
+ /* log bounds */
+ vl[0][0]= -3;vl[0][1]=5+1e-5;vl[0][2]=.1; /* 0..100 layer log bounds and spacing */
+ vl[1][0]= -3;vl[1][1]=5+1e-5;vl[1][2]=.1; /* 100..410 */
+ vl[2][0]= -3;vl[2][1]=5+1e-5;vl[2][2]=.1; /* 660 ... 2871 */
+ /* linear, in units of reference */
+ p->elayer[1] = 1.0; /* 410..660 reference, remains unchanged */
+ /* */
+ solved=0;
+ for(v1=vl[0][0];v1 <= vl[0][1];v1 += vl[0][2])
+ for(v2=vl[1][0];v2 <= vl[1][1];v2 += vl[1][2])
+ for(v3=vl[2][0];v3 <= vl[2][1];v3 += vl[2][2]){
+ /* layer viscosity structure */
+ p->elayer[0] = pow(10,v3);p->elayer[2] = pow(10,v2);p->elayer[3]=pow(10,v1);
+ hc_assign_viscosity(model,HC_INIT_E_FOUR_LAYERS,p->elayer,p);
+ /* print viscosities of 0...100, 100...410, and 660...2871 layer */
+ fprintf(stdout,"%14.7e %14.7e %14.7e\t",p->elayer[3],p->elayer[2],p->elayer[0]);
+ /* compute solution */
+ hc_solve(model,p->free_slip,p->solution_mode,sol_spectral,
+ (solved)?(FALSE):(TRUE), /* density changed? */
+ (solved)?(FALSE):(TRUE), /* plate velocity changed? */
+ TRUE, /* viscosity changed */
+ FALSE,p->compute_geoid,geoid,
+ p->verbose);
+ /* only output are the geoid correlations, for now */
+ hc_compute_correlation(geoid,p->ref_geoid,corr,1,p->verbose);
+ fprintf(stdout,"%10.7f %10.7f ",corr[0],corr[1]);
+ fprintf(stdout,"\n");
+ solved++;
+ }
+ break;
+ default:
+ HC_ERROR("hc","solver mode undefined");
}
+
/*
free memory
Modified: mc/1D/hc/trunk/make_tar
===================================================================
--- mc/1D/hc/trunk/make_tar 2008-05-24 02:57:03 UTC (rev 12019)
+++ mc/1D/hc/trunk/make_tar 2008-05-25 02:03:20 UTC (rev 12020)
@@ -2,7 +2,7 @@
#
# create a tar file
#
-ver=${1-0.1.4}
+ver=${1-0.1.5}
date=`date '+%m%d%y'`
tf=$HOME/tmp/hc-$ver.$date.tgz
Modified: mc/1D/hc/trunk/rick_sh_c.c
===================================================================
--- mc/1D/hc/trunk/rick_sh_c.c 2008-05-24 02:57:03 UTC (rev 12019)
+++ mc/1D/hc/trunk/rick_sh_c.c 2008-05-25 02:03:20 UTC (rev 12020)
@@ -464,8 +464,8 @@
// Legendre functions are precomputed
// */
double dpdt,dpdp,mphi,sin_theta,sfac,cfac;
- float *loc_plma=NULL,*loc_plmb=NULL,*plm=NULL,*dplm=NULL;
- int i,j,k,k2,m,l,lmaxp1,lm1;
+ float *plm=NULL,*dplm=NULL;
+ int i,k,k2,m,l,lmaxp1,lm1;
if(!rick->initialized){
fprintf(stderr,"rick_shc2d_pre_reg: error: initialize modules first\n");
exit(-1);
Modified: mc/1D/hc/trunk/sh_exp.c
===================================================================
--- mc/1D/hc/trunk/sh_exp.c 2008-05-24 02:57:03 UTC (rev 12019)
+++ mc/1D/hc/trunk/sh_exp.c 2008-05-25 02:03:20 UTC (rev 12020)
@@ -266,7 +266,6 @@
/* compute total correlation up to llim */
HC_PREC sh_correlation(struct sh_lms *exp1, struct sh_lms *exp2, int llim)
{
- HC_PREC corr;
return sh_correlation_per_degree(exp1,exp2,1,llim);
}
More information about the cig-commits
mailing list