[cig-commits] r20526 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Mon Jul 16 05:18:59 PDT 2012
Author: becker
Date: 2012-07-16 05:18:58 -0700 (Mon, 16 Jul 2012)
New Revision: 20526
Modified:
mc/1D/hc/trunk/Makefile
mc/1D/hc/trunk/Makefile.include
mc/1D/hc/trunk/ggrd_test.c
mc/1D/hc/trunk/hc.h
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_extract_spatial.c
mc/1D/hc/trunk/hc_init.c
mc/1D/hc/trunk/hc_matrix.c
mc/1D/hc/trunk/hc_misc.c
mc/1D/hc/trunk/hc_output.c
mc/1D/hc/trunk/hc_polsol.c
mc/1D/hc/trunk/hc_propagator.c
mc/1D/hc/trunk/hc_solve.c
mc/1D/hc/trunk/hc_torsol.c
mc/1D/hc/trunk/main.c
mc/1D/hc/trunk/prem_util.c
mc/1D/hc/trunk/rick_fft_c.c
mc/1D/hc/trunk/rick_sh_c.c
mc/1D/hc/trunk/sh.h
mc/1D/hc/trunk/sh_ana.c
mc/1D/hc/trunk/sh_corr.c
mc/1D/hc/trunk/sh_exp.c
mc/1D/hc/trunk/sh_model.c
mc/1D/hc/trunk/sh_power.c
mc/1D/hc/trunk/sh_rick.h
mc/1D/hc/trunk/sh_syn.c
mc/1D/hc/trunk/simple_test.c
Log:
Allow for quadruple precision.
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/Makefile 2012-07-16 12:18:58 UTC (rev 20526)
@@ -122,7 +122,7 @@
-L$(ODIR)/
#
-INC_FLAGS = $(HEAL_INC_FLAGS) \
+INC_FLAGS = $(HEAL_INC_FLAGS) $(ADD_FLAGS) \
$(RICK_INC_FLAGS) $(GGRD_INC_FLAGS)
#
# includes
@@ -258,7 +258,7 @@
mkdir -p $(BDIR);
clean:
- rm -f $(ODIR)/*.o $(ODIR)/*.a
+ rm -f $(ODIR)/*.o $(ODIR)/*.a $(BDIR)/*
#
# library
Modified: mc/1D/hc/trunk/Makefile.include
===================================================================
--- mc/1D/hc/trunk/Makefile.include 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/Makefile.include 2012-07-16 12:18:58 UTC (rev 20526)
@@ -6,5 +6,6 @@
#
# for GMT version >= 4.1.2, uncomment the next two lines
#
+#ADD_FLAGS = -Werror
GGRD_INC_FLAGS = -I$(GMTHOME)/include -I$(NETCDFHOME)/include
GGRD_LIBS_LINKLINE = -lggrd -lgmt -lpsl -lnetcdf
Modified: mc/1D/hc/trunk/ggrd_test.c
===================================================================
--- mc/1D/hc/trunk/ggrd_test.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/ggrd_test.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -3,7 +3,8 @@
int main(int argc, char **argv)
{
struct ggrd_master *ggrd;
- HC_PREC xloc[3],time,vr[4],vphi[4],vtheta[4],dtrange;
+ double time,vr[4],vphi[4],vtheta[4],dtrange;
+ double xloc[3];
static int order = 3;
hc_boolean calc_derivatives ;
double lon,lat,age;
@@ -54,7 +55,8 @@
/*
interpolate
*/
- if(ggrd_find_vel_and_der(xloc,time,dtrange,ggrd,order,calc_derivatives,
+ if(ggrd_find_vel_and_der(xloc,time,dtrange,ggrd,
+ order,calc_derivatives,
TRUE,vr,vtheta,vphi))
exit(-1);
if(interpolate_seafloor_ages(xloc[HC_THETA],
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/hc.h 2012-07-16 12:18:58 UTC (rev 20526)
@@ -26,24 +26,62 @@
#ifndef hc_boolean
#define hc_boolean unsigned short
#endif
-#ifndef HC_PREC /*
- precision for most C functions
- */
+
+#ifndef HC_PRECISION /* not defined */
+/*
+ default precision for most C functions
+*/
+#define HC_PRECISION 16 /* 8,16, or 32 denoting float, double,
+ or quad only for the purposes of
+ switching */
+#endif
+
+#ifndef __HC_PREC_DEFINED__
+
+#if HC_PRECISION == 32 /* quad prec */
+
+#define HC_PREC long double
+#define HC_HIGH_PREC long double
+#define HC_FLT_FORMAT "%Lf"
+#define HC_TWO_FLT_FORMAT "%Lf %Lf"
+#define HC_THREE_FLT_FORMAT "%Lf %Lf %Lf"
+#define HC_EPS_PREC 3e-25
+
+#elif HC_PRECISION == 16 /* double prec */
+
#define HC_PREC double
+#define HC_HIGH_PREC double /* some computations at higher prec */
#define HC_FLT_FORMAT "%lf"
#define HC_TWO_FLT_FORMAT "%lf %lf"
-#define HC_EPS_PREC 5e-15
+#define HC_THREE_FLT_FORMAT "%lf %lf %lf"
+#define HC_EPS_PREC 3e-16
+
+#else /* single prec */
+
+#define HC_PREC float
+#define HC_HIGH_PREC double
+#define HC_FLT_FORMAT "%f"
+#define HC_TWO_FLT_FORMAT "%f %f"
+#define HC_THREE_FLT_FORMAT "%f %f %f"
+#define HC_EPS_PREC 3e-7
+
#endif
+#define HC_BIN_PREC float /* precision for binary I/O */
+
#ifndef HC_CPREC
#define HC_CPREC HC_PREC
#endif
+#define __HC_PREC_DEFINED__
+#endif
+
+
#define HC_CHAR_LENGTH 300 /* length of char arrays */
-#define HC_BIN_PREC float /* precision for binary I/O */
+
#ifndef __HC_DEF_COMPLEX__
/* complex variables */
struct hc_scmplx{ /* single precision */
@@ -130,8 +168,8 @@
struct hc_ps{
/* scaling factors will only be computed once */
int ncalled;
- double rho_scale;
- double alpha, beta, geoid_factor;
+ HC_HIGH_PREC rho_scale;
+ HC_HIGH_PREC alpha, beta, geoid_factor;
hc_boolean rho_init,
prop_params_init, /* parameters for propagator computation */
abg_init , /* alpha, beta factors */
@@ -288,10 +326,10 @@
*/
/* poloidal solution */
struct sh_lms *pol_sol;
- double *rho,*rho_zero; /*
+ HC_PREC *rho,*rho_zero; /*
density factors
*/
- double *rprops,*pvisc,*props,*ppots, /* propagator related */
+ HC_HIGH_PREC *rprops,*pvisc,*props,*ppots, /* propagator related */
*den;
/*
propagator related factors as well
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2012-07-16 12:18:58 UTC (rev 20526)
@@ -60,6 +60,7 @@
void hc_ludcmp_3x3(double [3][3], int, int *);
void hc_lubksb_3x3(double [3][3], int, int *, double *);
/* hc_misc.c */
+void hc_hvecalloc(double **, int, char *);
void hc_dvecalloc(double **, int, char *);
void hc_svecalloc(float **, int, char *);
void hc_ivecalloc(int **, int, char *);
@@ -68,8 +69,8 @@
void hc_svecrealloc(float **, int, char *);
void hc_dvecrealloc(double **, int, char *);
void hc_vecrealloc(double **, int, char *);
-float hc_svec_rms_diff(float *, float *, int);
-float hc_svec_rms(float *, int);
+float hc_vec_rms_diff(double *, double *, int);
+float hc_vec_rms(double *, int);
void hc_a_equals_b_svector(float *, float *, int);
void hc_a_equals_b_vector(double *, double *, int);
float hc_mean_svec(float *, int);
@@ -81,14 +82,14 @@
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);
-void lonlatpv2cv(float, float, float *, float *);
-void lonlatpv2cv_with_base(float *, double *, float *);
-void calc_polar_base_at_theta_phi(float, float, double *);
+void lonlatpv2cv(double, float, double *, double *);
+void lonlatpv2cv_with_base(double *, double *, double *);
+void calc_polar_base_at_theta_phi(double, double, double *);
void hc_linear_interpolate(double *, int, double, int *, int *, double *, double *);
/* 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);
-void hc_print_spatial_solution(struct hcs *, struct sh_lms *, float *, char *, char *, int, unsigned short, unsigned short);
+void hc_print_spatial_solution(struct hcs *, struct sh_lms *, double *, char *, char *, int, unsigned short, unsigned short);
void hc_print_depth_layers(struct hcs *, FILE *, unsigned short);
void hc_print_3x3(double [3][3], FILE *);
void hc_print_sm(double [6][4], FILE *);
@@ -99,8 +100,10 @@
void hc_compute_solution_scaling_factors(struct hcs *, int, double, double, double *);
void hc_print_poloidal_solution(struct sh_lms *, struct hcs *, int, char *, unsigned short, unsigned short);
void hc_print_toroidal_solution(double *, int, struct hcs *, int, char *, unsigned short);
-void hc_print_vtk(FILE *, float *, float *, int, int, unsigned short, int, float *, int, int);
-void hc_print_be_float(float *, int, FILE *, unsigned short);
+void hc_print_vtk(FILE *, double *, double *, int, int, unsigned short, int, double *, int, int);
+int hc_print_be_float(double *, int, FILE *, unsigned short);
+int hc_print_float(double *, int, FILE *);
+int hc_read_float(double *, int, FILE *);
void hc_print_be_int(int *, int, FILE *, unsigned short);
unsigned short hc_is_little_endian(void);
void hc_flip_byte_order(void *, size_t);
@@ -114,7 +117,7 @@
/* hc_solve.c */
void hc_solve(struct hcs *, unsigned short, int, struct sh_lms *, unsigned short, unsigned short, unsigned short, unsigned short, unsigned short, struct sh_lms *, struct sh_lms *, struct sh_lms *, unsigned short);
void hc_sum(struct hcs *, int, struct sh_lms *, struct sh_lms *, int, unsigned short, struct sh_lms *, unsigned short);
-void hc_compute_sol_spatial(struct hcs *, struct sh_lms *, float **, unsigned short);
+void hc_compute_sol_spatial(struct hcs *, struct sh_lms *, double **, unsigned short);
/* hc_torsol.c */
void hc_torsol(struct hcs *, int, int, int, double *, double **, double **, struct sh_lms *, struct sh_lms *, double *, unsigned short);
/* main.c */
@@ -132,25 +135,25 @@
int prem_read_para_set(double *, int, int, FILE *);
/* print_gauss_lat.c */
/* rick_fft_c.c */
-void rick_cs2ab(float *, int);
-void rick_ab2cs(float *, int);
-void rick_realft_nr(float *, int, int);
-void rick_four1_nr(float *, int, int);
+void rick_cs2ab(double *, int);
+void rick_ab2cs(double *, int);
+void rick_realft_nr(double *, int, int);
+void rick_four1_nr(double *, int, int);
/* rick_sh_c.c */
-void rick_compute_allplm(int, int, float *, float *, struct rick_module *);
-void rick_compute_allplm_reg(int, int, float *, float *, struct rick_module *, float *, int);
-void rick_pix2ang(int, int, float *, float *, struct rick_module *);
-void rick_shc2d(float *, float *, int, int, float *, float *, struct rick_module *);
-void rick_shc2d_reg(float *, float *, int, int, float *, float *, struct rick_module *, float *, int, float *, int, unsigned short);
-void rick_shc2d_pre(float *, float *, int, float *, float *, int, float *, float *, struct rick_module *);
-void rick_shc2d_pre_reg(float *, float *, int, float *, float *, int, float *, float *, struct rick_module *, float *, int, float *, int, unsigned short);
-void rick_shc2d_irreg(float *, float *, int, int, float *, float *, struct rick_module *, float *, float *, int);
-void rick_shd2c(float *, float *, int, int, float *, float *, struct rick_module *);
-void rick_shd2c_pre(float *, float *, int, float *, float *, int, float *, float *, struct rick_module *);
+void rick_compute_allplm(int, int, double *, double *, struct rick_module *);
+void rick_compute_allplm_reg(int, int, double *, double *, struct rick_module *, double *, int);
+void rick_pix2ang(int, int, double *, double *, struct rick_module *);
+void rick_shc2d(double *, double *, int, int, double *, double *, struct rick_module *);
+void rick_shc2d_reg(double *, double *, int, int, double *, double *, struct rick_module *, double *, int, double *, int, unsigned short);
+void rick_shc2d_pre(double *, double *, int, double *, double *, int, double *, double *, struct rick_module *);
+void rick_shc2d_pre_reg(double *, double *, int, double *, double *, int, double *, double *, struct rick_module *, double *, int, double *, int, unsigned short);
+void rick_shc2d_irreg(double *, double *, int, int, double *, double *, struct rick_module *, double *, double *, int);
+void rick_shd2c(double *, double *, int, int, double *, double *, struct rick_module *);
+void rick_shd2c_pre(double *, double *, int, double *, double *, int, double *, double *, struct rick_module *);
void rick_init(int, int, int *, int *, int *, struct rick_module *, unsigned short);
void rick_free_module(struct rick_module *, int);
-void rick_plmbar1(float *, float *, int, int, float, struct rick_module *);
-void rick_gauleg(float, float, float *, float *, int);
+void rick_plmbar1(double *, double *, int, int, double, struct rick_module *);
+void rick_gauleg(double, double, double *, double *, int);
/* rotvec2vel.c */
FILE *rv_myopen(const char *, const char *);
/* sh_ana.c */
@@ -171,7 +174,7 @@
void sh_free_expansion(struct sh_lms *, int);
void sh_clear_alm(struct sh_lms *);
double sh_total_power(struct sh_lms *);
-void sh_compute_power_per_degree(struct sh_lms *, float *);
+void sh_compute_power_per_degree(struct sh_lms *, double *);
double sh_correlation(struct sh_lms *, struct sh_lms *, int);
double sh_correlation_per_degree(struct sh_lms *, struct sh_lms *, int, int);
void sh_print_parameters_to_file(struct sh_lms *, int, int, int, double, FILE *, unsigned short, unsigned short, unsigned short);
@@ -179,22 +182,22 @@
void sh_print_coefficients_to_file(struct sh_lms *, int, FILE *, double *, unsigned short, unsigned short);
void sh_read_coefficients_from_file(struct sh_lms *, int, int, FILE *, unsigned short, double *, unsigned short);
void sh_print_nonzero_coeff(struct sh_lms *, FILE *);
-void sh_read_spatial_data_from_file(struct sh_lms *, FILE *, unsigned short, int, float *, float *);
-void sh_read_spatial_data_from_grd(struct sh_lms *, struct ggrd_gt *, unsigned short, int, float *, float *);
-void sh_read_spatial_data(struct sh_lms *, FILE *, struct ggrd_gt *, unsigned short, unsigned short, int, float *, float *);
-void sh_compute_spatial_basis(struct sh_lms *, FILE *, unsigned short, float, float **, int, unsigned short);
-void sh_compute_spectral(float *, int, unsigned short, float **, struct sh_lms *, unsigned short);
-void sh_compute_spatial(struct sh_lms *, int, unsigned short, float **, float *, unsigned short);
-void sh_compute_spatial_reg(struct sh_lms *, int, unsigned short, float **, float *, int, float *, int, float *, unsigned short, unsigned short);
-void sh_compute_spatial_irreg(struct sh_lms *, int, float *, float *, int, float *, unsigned short);
+void sh_read_spatial_data_from_file(struct sh_lms *, FILE *, unsigned short, int, double *, double *);
+void sh_read_spatial_data_from_grd(struct sh_lms *, struct ggrd_gt *, unsigned short, int, double *, double *);
+void sh_read_spatial_data(struct sh_lms *, FILE *, struct ggrd_gt *, unsigned short, unsigned short, int, double *, double *);
+void sh_compute_spatial_basis(struct sh_lms *, FILE *, unsigned short, double, double **, int, unsigned short);
+void sh_compute_spectral(double *, int, unsigned short, double **, struct sh_lms *, unsigned short);
+void sh_compute_spatial(struct sh_lms *, int, unsigned short, double **, double *, unsigned short);
+void sh_compute_spatial_reg(struct sh_lms *, int, unsigned short, double **, double *, int, double *, int, double *, unsigned short, unsigned short);
+void sh_compute_spatial_irreg(struct sh_lms *, int, double *, double *, int, double *, unsigned short);
void sh_exp_type_error(char *, struct sh_lms *);
-void sh_print_plm(float *, int, int, int, FILE *);
-void sh_print_spatial_data_to_file(struct sh_lms *, int, float *, unsigned short, float, FILE *);
-void sh_get_coordinates(struct sh_lms *, int, float *, float *);
-void sh_print_reg_spatial_data_to_file(struct sh_lms *, int, float *, unsigned short, float, float *, int, float *, int, FILE *);
-void sh_print_irreg_spatial_data_to_file(struct sh_lms *, int, float *, unsigned short, float, float *, float *, int, FILE *);
-void sh_compute_plm(struct sh_lms *, int, float **, unsigned short);
-void sh_compute_plm_reg(struct sh_lms *, int, float **, unsigned short, float *, int);
+void sh_print_plm(double *, int, int, int, FILE *);
+void sh_print_spatial_data_to_file(struct sh_lms *, int, double *, unsigned short, double, FILE *);
+void sh_get_coordinates(struct sh_lms *, int, double *, double *);
+void sh_print_reg_spatial_data_to_file(struct sh_lms *, int, double *, unsigned short, double, double *, int, double *, int, FILE *);
+void sh_print_irreg_spatial_data_to_file(struct sh_lms *, int, double *, unsigned short, double, double *, double *, int, FILE *);
+void sh_compute_plm(struct sh_lms *, int, double **, unsigned short);
+void sh_compute_plm_reg(struct sh_lms *, int, double **, unsigned short, double *, int);
void sh_get_coeff(struct sh_lms *, int, int, int, unsigned short, double *);
void sh_write_coeff(struct sh_lms *, int, int, int, unsigned short, double *);
void sh_add_coeff(struct sh_lms *, int, int, int, unsigned short, double *);
@@ -209,14 +212,13 @@
void sh_free_model(struct sh_lms_model *);
void sh_print_model_coefficients(struct sh_lms_model *, FILE *, unsigned short, unsigned short);
void sh_print_model_spatial_basis(struct sh_lms_model *, FILE *, unsigned short);
-void sh_read_model_spatial_data(struct sh_lms_model *, float **, FILE *, unsigned short);
-void sh_compute_model_spectral(struct sh_lms_model *, float *, unsigned short);
-void sh_compute_model_spatial(struct sh_lms_model *, float **, unsigned short);
-void sh_print_model_spatial_data(struct sh_lms_model *, float *, FILE *, unsigned short);
+void sh_read_model_spatial_data(struct sh_lms_model *, double **, FILE *, unsigned short);
+void sh_compute_model_spectral(struct sh_lms_model *, double *, unsigned short);
+void sh_compute_model_spatial(struct sh_lms_model *, double **, unsigned short);
+void sh_print_model_spatial_data(struct sh_lms_model *, double *, FILE *, unsigned short);
/* sh_power.c */
/* sh_syn.c */
/* sh_test.c */
/* simple_test.c */
-void hc_dvecalloc(double **, int, char *);
/* spherepack_sh.c */
/* test_fft.c */
Modified: mc/1D/hc/trunk/hc_extract_spatial.c
===================================================================
--- mc/1D/hc/trunk/hc_extract_spatial.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/hc_extract_spatial.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -17,9 +17,9 @@
struct hcs *model;
HC_PREC zlabel;
hc_boolean binary_in = TRUE, verbose = FALSE,read_dsol=FALSE;
- float *data,*plm=NULL,*xpos,*xvec,lon,lat,theta,phi,xtmp[3],pvec[3],
+ HC_PREC *data,*plm=NULL,*xpos,*xvec,lon,lat,theta,phi,xtmp[3],pvec[3],
*xscalar;
- double polar_base[9];
+ HC_PREC polar_base[9];
hc_struc_init(&model);
/*
deal with parameters
@@ -132,9 +132,9 @@
ndata_all = npoints * (shps + shps_read_d);
if((mode == 5)||(mode==6)){ /* save all layers */
- hc_svecalloc(&data,model->nradp2 * ndata_all,"hc_extract_spatial");
+ hc_vecalloc(&data,model->nradp2 * ndata_all,"hc_extract_spatial");
}else
- hc_svecalloc(&data, ndata_all,"hc_extract_spatial");
+ hc_vecalloc(&data, ndata_all,"hc_extract_spatial");
for(lc=0,ilayer=i1;ilayer <= i2;ilayer++,lc++){
/*
output
@@ -195,10 +195,10 @@
*/
if(shps != 3)HC_ERROR("hc_extract_spatial","shps has to be 3 for mode 5 and 6");
/* convert */
- hc_svecalloc(&xpos,model->nradp2 * ndata,"hc_extract_spatial");
- hc_svecalloc(&xvec,model->nradp2 * ndata,"hc_extract_spatial");
+ hc_vecalloc(&xpos,model->nradp2 * ndata,"hc_extract_spatial");
+ hc_vecalloc(&xvec,model->nradp2 * ndata,"hc_extract_spatial");
if(read_dsol)
- hc_svecalloc(&xscalar,model->nradp2 * ndata_d,"hc_extract_spatial");
+ hc_vecalloc(&xscalar,model->nradp2 * ndata_d,"hc_extract_spatial");
for(i=0;i < npoints;i++){ /* loop through all points */
/* lon lat coordinates */
sh_get_coordinates((vsol+i1*3),i,&lon,&lat);
Modified: mc/1D/hc/trunk/hc_init.c
===================================================================
--- mc/1D/hc/trunk/hc_init.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/hc_init.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -407,7 +407,7 @@
fprintf(stderr,"-vshs\t\tuse the short format (only lmax in header) for the plate velocities (%s)\n",
hc_name_boolean(p->read_short_pvel_sh));
fprintf(stderr,"-vdir\t\tvelocities are given in files name/vel.1.ab to vel.%i.ab for different times,\n\t\t-%g to -1 Ma before present, where name is from -pvel\n",
- HC_PVEL_TSTEPS,(double)HC_PVEL_TSTEPS);
+ HC_PVEL_TSTEPS,(HC_PREC)HC_PVEL_TSTEPS);
fprintf(stderr,"-vtime\ttime\tuse this particular time step of the plate velocities (%g)\n\n",
p->pvel_time);
@@ -581,7 +581,7 @@
hc_get_flt_frmt_string(fstring,2,FALSE);
rold = hc->r_cmb;
/* start read loop */
- while(fscanf(in,"%lf %lf",(hc->rvisc+hc->nvis),(hc->visc+hc->nvis))==2){
+ while(fscanf(in,HC_TWO_FLT_FORMAT,(hc->rvisc+hc->nvis),(hc->visc+hc->nvis))==2){
if(hc->visc[hc->nvis] < 1e15)
fprintf(stderr,"hc_assign_viscosity: WARNING: expecting viscosities in Pas, read %g at layer %i\n",
hc->visc[hc->nvis],hc->nvis);
@@ -678,9 +678,10 @@
{
FILE *in;
int type,lmax,shps,ilayer,nset,ivec,i,j;
- HC_PREC *dtop,*dbot,zlabel,local_scale,dens_scale[1],rho0;
+ HC_PREC *dtop,*dbot,zlabel,local_scale,dens_scale[1];
+ double rho0;
hc_boolean reported = FALSE,read_on;
- double dtmp[3];
+ HC_PREC dtmp[3];
hc->compressible = compressible;
hc->inho = 0;
if(hc->dens_init) /* clear old expansions, if
@@ -737,7 +738,7 @@
while(read_on){
if(use_short_format){
/* short format I/O */
- i = fscanf(in,"%lf",dtmp);zlabel = (HC_PREC)dtmp[0];
+ i = fscanf(in,HC_FLT_FORMAT,dtmp);zlabel = (HC_PREC)dtmp[0];
i += fscanf(in,"%i",&lmax);
read_on = (i == 2)?(TRUE):(FALSE);
ivec = 0;shps = 1;type = HC_DEFAULT_INTERNAL_FORMAT;
@@ -777,7 +778,7 @@
/*
assign depth, this assumes that we are reading in depths [km]
*/
- hc->rden[hc->inho] = HC_ND_RADIUS((double)zlabel);
+ hc->rden[hc->inho] = HC_ND_RADIUS((HC_PREC)zlabel);
if(scale_dens_anom_with_prem){
/*
@@ -786,25 +787,26 @@
*/
- prem_get_rho(&rho0,hc->rden[hc->inho],hc->prem);
+ prem_get_rho(&rho0,(double)(hc->rden[hc->inho]),
+ hc->prem);
+
rho0 /= 1000.0;
if(rho0 < 3)
fprintf(stderr,"\nhc_assign_density: WARNING: using small (%g) density from PREM for layer at depth %g\n\n",
rho0*1000,HC_Z_DEPTH(hc->rden[hc->inho]));
}else{
/* mean value */
- rho0 = hc->avg_den_mantle;
+ rho0 = (double)hc->avg_den_mantle;
}
/*
density anomaly
*/
/* scaling factor without depth dependence */
- dens_scale[0] = ((HC_PREC)HC_DENSITY_SCALING ) * rho0;
+ dens_scale[0] = HC_DENSITY_SCALING * (HC_PREC)rho0;
if(verbose >= 2){
-
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/ hc->avg_den_mantle,dens_scale[0],hc->inho+1,nset,zlabel);
+ (double)hc->rden[hc->inho],
+ HC_DENSITY_SCALING,rho0/ (double)hc->avg_den_mantle,(double)dens_scale[0],hc->inho+1,nset,(double)zlabel);
}
if(hc->inho){
/*
@@ -879,7 +881,7 @@
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],
+ (double)hc->rden[i],
(dd_dens_scale == HC_DD_READ_FROM_FILE)?("depth-dependent"):((dd_dens_scale==HC_DD_CONSTANT)?("constant"):("polynomial")),local_scale,i,hc->inho);
}
}
@@ -900,7 +902,7 @@
hc->r[0] = hc->r_cmb; /* CMB */
if(hc->rden[0] <= hc->r[0]){
- fprintf(stderr,"hc_assign_density: rden[0]: %g r[0]: %g\n",hc->rden[0], hc->r[0]);
+ fprintf(stderr,"hc_assign_density: rden[0]: %g r[0]: %g\n",(double)hc->rden[0],(double) hc->r[0]);
HC_ERROR("hc_assign_density","first density layer has to be above internal CMB limit");
}
for(i=0;i<hc->nrad;i++) /* density layers */
@@ -950,8 +952,8 @@
if(verbose)
for(i=0;i < hc->nrad;i++)
fprintf(stderr,"hc_assign_density: dens %3i: r: %8.6f df: %8.6f |rho|: %8.4f\n",
- i+1,hc->rden[i],hc->dfact[i],
- sqrt(sh_total_power((hc->dens_anom+i))));
+ i+1,(double)hc->rden[i],(double)hc->dfact[i],
+ (double)sqrt(sh_total_power((hc->dens_anom+i))));
free(dbot);free(dtop);
} /* end layer structure part */
@@ -1248,7 +1250,7 @@
HC_PREC rcmb)
{
HC_PREC smean;
- double dtmp[2];
+ HC_PREC dtmp[2];
int i;
FILE *in;
if(p->dd_dens_scale == HC_DD_READ_FROM_FILE){
@@ -1266,7 +1268,7 @@
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){
+ while(fscanf(in,HC_TWO_FLT_FORMAT,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];
@@ -1307,7 +1309,7 @@
void hc_read_geoid(struct hc_parameters *p)
{
int type,lmax,shps,ilayer,nset,ivec;
- double zlabel;
+ HC_PREC zlabel;
FILE *in;
HC_PREC fac[3]={1,1,1};
Modified: mc/1D/hc/trunk/hc_matrix.c
===================================================================
--- mc/1D/hc/trunk/hc_matrix.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/hc_matrix.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -34,10 +34,10 @@
for(i=0;i<n;i++)
for(j=0;j<n;j++)
- fscanf(stdin,"%lf",&amat[i][j]);
+ fscanf(stdin,HC_FLT_FORMAT,&amat[i][j]);
/* read in x from stdin */
for(i=0;i<n;i++)
- fscanf(stdin,"%lf",&bvec[i]);
+ fscanf(stdin,HC_FLT_FORMAT,&bvec[i]);
/* solve A x = b, where b will be modified */
hc_ludcmp_3x3(amat,n,indx);hc_lubksb_3x3(amat,n,indx,bvec);
for(i=0;i<n;i++)fprintf(stdout,"%g\n",bvec[i]);
Modified: mc/1D/hc/trunk/hc_misc.c
===================================================================
--- mc/1D/hc/trunk/hc_misc.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/hc_misc.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -12,7 +12,14 @@
-/* double vector allocation */
+/* high precision vector allocation */
+void hc_hvecalloc(HC_HIGH_PREC **x,int n,char *message)
+{
+ *x = (HC_HIGH_PREC *)malloc(sizeof(HC_HIGH_PREC)*(size_t)n);
+ if(! (*x))
+ HC_MEMERROR(message);
+}
+/* double */
void hc_dvecalloc(double **x,int n,char *message)
{
*x = (double *)malloc(sizeof(double)*(size_t)n);
@@ -26,7 +33,7 @@
if(! (*x))
HC_MEMERROR(message);
}
-/* integer vector allocation */
+/* integer vector allocation */
void hc_ivecalloc(int **x,int n,char *message)
{
*x = (int *)malloc(sizeof(int)*(size_t)n);
@@ -56,10 +63,10 @@
if(!(*x))
HC_MEMERROR(message);
}
-/* double vector reallocation */
-void hc_dvecrealloc(double **x,int n,char *message)
+/* HC_HIGH_PREC vector reallocation */
+void hc_dvecrealloc(HC_HIGH_PREC **x,int n,char *message)
{
- *x = (double *)realloc(*x,sizeof(double)*(size_t)n);
+ *x = (HC_HIGH_PREC *)realloc(*x,sizeof(HC_HIGH_PREC)*(size_t)n);
if(!(*x))
HC_MEMERROR(message);
}
@@ -73,20 +80,20 @@
/*
sqrt(sum(squared diff)) between [n] vectors a and b
*/
-float hc_svec_rms_diff(float *a,float *b,int n)
+float hc_vec_rms_diff(HC_PREC *a,HC_PREC *b,int n)
{
int i;
- float sum=0.0,tmp;
+ HC_PREC sum=0.0,tmp;
for(i=0;i<n;i++){
tmp = a[i] - b[i];
sum += tmp*tmp;
}
return sqrt(sum/n);
}
-float hc_svec_rms(float *a,int n)
+float hc_vec_rms(HC_PREC *a,int n)
{
int i;
- float sum=0.0;
+ HC_PREC sum=0.0;
for(i=0;i<n;i++){
sum += a[i] * a[i];
}
@@ -125,8 +132,8 @@
return sum;
}
-/* zero a double precision vector */
-void hc_zero_dvector(double *x, int n)
+/* zero a HC_HIGH_PREC precision vector */
+void hc_zero_dvector(HC_HIGH_PREC *x, int n)
{
int i;
for(i=0;i<n;i++)
@@ -159,6 +166,8 @@
sprintf(type_s,"f");
}else if (sizeof(HC_PREC) == sizeof(double)){
sprintf(type_s,"lf");
+ }else if (sizeof(HC_PREC) == sizeof(long double)){
+ sprintf(type_s,"lf");
}else{
fprintf(stderr,"hc_get_flt_frmt_string: assignment error\n");
exit(-1);
@@ -233,19 +242,19 @@
*/
-void lonlatpv2cv(float lon, float lat,
- float *polar_vec,float *cart_vec)
+void lonlatpv2cv(HC_PREC lon, float lat,
+ HC_PREC *polar_vec,HC_PREC *cart_vec)
{
- double polar_base[9];
- float theta,phi;
- theta = LAT2THETA((double)lat);
- phi = LON2PHI((double)lon);
+ HC_HIGH_PREC polar_base[9];
+ HC_PREC theta,phi;
+ theta = LAT2THETA((HC_HIGH_PREC)lat);
+ phi = LON2PHI((HC_HIGH_PREC)lon);
calc_polar_base_at_theta_phi(theta,phi,polar_base);
lonlatpv2cv_with_base(polar_vec,polar_base,cart_vec);
}
-void lonlatpv2cv_with_base(float *polar_vec,
- double *polar_base,
- float *cart_vec)
+void lonlatpv2cv_with_base(HC_PREC *polar_vec,
+ HC_HIGH_PREC *polar_base,
+ HC_PREC *cart_vec)
{
int i;
// convert vector
@@ -255,13 +264,13 @@
cart_vec[i] += polar_base[6+i] * polar_vec[2];
}
}
-void calc_polar_base_at_theta_phi(float theta, float phi,
- double *polar_base)
+void calc_polar_base_at_theta_phi(HC_PREC theta, HC_PREC phi,
+ HC_HIGH_PREC *polar_base)
{
- double cp,sp,ct,st;
+ HC_HIGH_PREC cp,sp,ct,st;
// base vecs
- ct=cos((double)theta);cp=cos((double)phi);
- st=sin((double)theta);sp=sin((double)phi);
+ ct=cos((HC_HIGH_PREC)theta);cp=cos((HC_HIGH_PREC)phi);
+ st=sin((HC_HIGH_PREC)theta);sp=sin((HC_HIGH_PREC)phi);
//
polar_base[0]= st * cp;
polar_base[1]= st * sp;
Modified: mc/1D/hc/trunk/hc_output.c
===================================================================
--- mc/1D/hc/trunk/hc_output.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/hc_output.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -15,11 +15,13 @@
*/
-void hc_print_spectral_solution(struct hcs *hc,struct sh_lms *sol,FILE *out,int sol_mode,
- hc_boolean binary, hc_boolean verbose)
+void hc_print_spectral_solution(struct hcs *hc,struct sh_lms *sol,
+ FILE *out,int sol_mode,
+ hc_boolean binary,
+ hc_boolean verbose)
{
int i,os;
- static int ntype = 3; /* three sets of solutions, r/pol/tor */
+ const int ntype = 3; /* three sets of solutions, r/pol/tor */
HC_PREC fac[3];
if(!hc->spectral_solution_computed)
HC_ERROR("hc_print_spectral_solution","spectral solution not computed");
@@ -32,7 +34,8 @@
scale to cm/yr, or MPa for stress solutions
*/
- hc_compute_solution_scaling_factors(hc,sol_mode,hc->r[i],hc->dvisc[i],fac);
+ hc_compute_solution_scaling_factors(hc,sol_mode,
+ hc->r[i],hc->dvisc[i],fac);
/*
write parameters, convert radius to depth in [km]
*/
@@ -50,24 +53,27 @@
switch(sol_mode){
case HC_VEL:
fprintf(stderr,"hc_print_spectral_solution: z: %8.3f vel: |r|: %11.3e |pol|: %11.3e |tor|: %11.3e (scale: %g cm/yr)\n",
- HC_Z_DEPTH(hc->r[i]),sqrt(sh_total_power((sol+os))),
- sqrt(sh_total_power((sol+os+1))),
- sqrt(sh_total_power((sol+os+2))),
- fac[0]/11.1194926644559);
+ (double)HC_Z_DEPTH(hc->r[i]),
+ (double)sqrt(sh_total_power((sol+os))),
+ (double)sqrt(sh_total_power((sol+os+1))),
+ (double)sqrt(sh_total_power((sol+os+2))),
+ (double)(fac[0]/11.1194926644559));
break;
case HC_RTRACTIONS:
fprintf(stderr,"hc_print_spectral_solution: z: %8.3f rtrac: |r|: %11.3e |pol|: %11.3e |tor|: %11.3e (scale: %g MPa)\n",
- HC_Z_DEPTH(hc->r[i]),sqrt(sh_total_power((sol+os))),
- sqrt(sh_total_power((sol+os+1))),
- sqrt(sh_total_power((sol+os+2))),
- fac[0]/(0.553073278428428/hc->r[i]));
+ (double)HC_Z_DEPTH(hc->r[i]),
+ (double)sqrt(sh_total_power((sol+os))),
+ (double)sqrt(sh_total_power((sol+os+1))),
+ (double)sqrt(sh_total_power((sol+os+2))),
+ (double)(fac[0]/(0.553073278428428/hc->r[i])));
break;
case HC_HTRACTIONS:
fprintf(stderr,"hc_print_spectral_solution: z: %8.3f htrac: |r|: %11.3e |pol|: %11.3e |tor|: %11.3e (scale: %g MPa)\n",
- HC_Z_DEPTH(hc->r[i]),sqrt(sh_total_power((sol+os))),
- sqrt(sh_total_power((sol+os+1))),
- sqrt(sh_total_power((sol+os+2))),
- fac[0]/(0.553073278428428/hc->r[i]));
+ (double)HC_Z_DEPTH(hc->r[i]),
+ (double)sqrt(sh_total_power((sol+os))),
+ (double)sqrt(sh_total_power((sol+os+1))),
+ (double)sqrt(sh_total_power((sol+os+2))),
+ (double)(fac[0]/(0.553073278428428/hc->r[i])));
break;
default:
fprintf(stderr,"hc_print_spectral_solution: sol mode %i undefined\n",sol_mode);
@@ -87,11 +93,14 @@
*/
-void hc_print_sh_scalar_field(struct sh_lms *sh, FILE *out, hc_boolean short_format,
- hc_boolean binary, hc_boolean verbose)
+void hc_print_sh_scalar_field(struct sh_lms *sh, FILE *out,
+ hc_boolean short_format,
+ hc_boolean binary,
+ hc_boolean verbose)
{
HC_CPREC fac[1] = {1.0};
- sh_print_parameters_to_file(sh,1,0,1,0.0,out,short_format,binary,verbose); /* parameters in long format */
+ sh_print_parameters_to_file(sh,1,0,1,0.0,out,
+ short_format,binary,verbose); /* parameters in long format */
sh_print_coefficients_to_file(sh,1,out,fac,binary,verbose); /* coefficients */
}
@@ -110,14 +119,14 @@
*/
void hc_print_spatial_solution(struct hcs *hc,
struct sh_lms *sol,
- float *sol_x, char *name,
+ HC_PREC *sol_x, char *name,
char *dfilename,
int sol_mode, hc_boolean binary,
hc_boolean verbose)
{
int i,j,k,os[3],los,np,np2,np3;
FILE *file_dummy=NULL,*out,*dout;
- float flt_dummy=0,*xy=NULL,value[3];
+ HC_PREC flt_dummy=0,*xy=NULL,value[3];
HC_PREC fac[3];
char filename[300];
if(!hc->spatial_solution_computed)
@@ -152,7 +161,7 @@
hc->dvisc[i],fac);
/* write depth in [km] to dout file */
- fprintf(dout,"%g\n",HC_Z_DEPTH(hc->r[i]));
+ fprintf(dout,"%g\n",(double)HC_Z_DEPTH(hc->r[i]));
for(k=0;k < 3;k++) /* pointers */
os[k] = i * np3 + k*np;
/*
@@ -172,10 +181,10 @@
sprintf(filename,"%s.%i.bin",name,i+1);
out = ggrd_open(filename,"w","hc_print_spatial_solution");
for(j=los=0;j < np;j++,los+=2){ /* loop through all points in layer */
- fwrite((xy+los),sizeof(float),2,out);
+ hc_print_float((xy+los),2,out);
for(k=0;k<3;k++)
value[k] = sol_x[os[k]] * fac[k];
- fwrite(value,sizeof(float),3,out);
+ hc_print_float(value,3,out);
os[0]++;os[1]++;os[2]++;
}
fclose(out);
@@ -187,15 +196,19 @@
for(k=0;k<3;k++)
value[k] = sol_x[os[k]] * fac[k];
fprintf(out,"%11g %11g\t%12.5e %12.5e %12.5e\n",
- xy[los],xy[los+1],value[0],value[1],value[2]);
+ (double)xy[los],
+ (double)xy[los+1],
+ (double)value[0],
+ (double)value[1],(double)value[2]);
os[0]++;os[1]++;os[2]++;
}
fclose(out);
}
if(verbose >= 2)
fprintf(stderr,"hc_print_spatial_solution: layer %3i: RMS: r: %12.5e t: %12.5e p: %12.5e file: %s\n",
- i+1,hc_svec_rms((sol_x+i*np3),np),hc_svec_rms((sol_x+i*np3+np),np),
- hc_svec_rms((sol_x+i*np3+np2),np),
+ i+1, (double)hc_vec_rms((sol_x+i*np3),np),
+ (double)hc_vec_rms((sol_x+i*np3+np),np),
+ (double)hc_vec_rms((sol_x+i*np3+np2),np),
filename);
}
fclose(dout);
@@ -216,7 +229,7 @@
int i;
/* number of solution sets of ntype solutions */
for(i=0;i < hc->nradp2;i++)
- fprintf(out,"%g\n",HC_Z_DEPTH(hc->r[i]));
+ fprintf(out,"%g\n",(double)HC_Z_DEPTH(hc->r[i]));
}
@@ -230,7 +243,7 @@
int i,j;
for(i=0;i<3;i++){
for(j=0;j<3;j++)
- fprintf(out,"%11.4e ",a[i][j]);
+ fprintf(out,"%11.4e ",(double)a[i][j]);
fprintf(out,"\n");
}
}
@@ -244,7 +257,7 @@
int i,j;
for(i=0;i < 6;i++){
for(j=0;j<4;j++)
- fprintf(out,"%11.4e ",a[i][j]);
+ fprintf(out,"%11.4e ",(double)a[i][j]);
fprintf(out,"\n");
}
}
@@ -253,7 +266,7 @@
{
int i;
for(i=0;i<n;i++)
- fprintf(out,"%11.4e ",a[i]);
+ fprintf(out,"%11.4e ",(double)a[i]);
fprintf(out,"\n");
}
void hc_print_vector_label(HC_PREC *a, int n,FILE *out,
@@ -262,7 +275,7 @@
int i;
fprintf(out,"%s: ",label);
for(i=0;i<n;i++)
- fprintf(out,"%11.4e ",a[i]);
+ fprintf(out,"%11.4e ",(double)a[i]);
fprintf(out,"\n");
}
void hc_print_matrix_label(HC_PREC *a, int m,
@@ -272,7 +285,7 @@
for(j=0;j<m;j++){
fprintf(out,"%s: ",label);
for(i=0;i<n;i++)
- fprintf(out,"%11.4e ",a[j*n+i]);
+ fprintf(out,"%11.4e ",(double)a[j*n+i]);
fprintf(out,"\n");
}
}
@@ -282,7 +295,7 @@
{
int i;
for(i=0;i<n;i++)
- fprintf(out,"%11.4e\n",a[i]);
+ fprintf(out,"%11.4e\n",(double)a[i]);
}
/*
@@ -345,10 +358,12 @@
alim = (m==0)?(1):(2);
for(a_or_b=0;a_or_b < alim;a_or_b++){
for(i=os=0;i < nl;i++,os+=6){
- fprintf(out,"%3i %3i %1i %3i %8.5f ",l,m,a_or_b,i+1,hc->r[i]);
+ fprintf(out,"%3i %3i %1i %3i %8.5f ",
+ l,m,a_or_b,i+1,(double)hc->r[i]);
for(j=0;j < 6;j++){
- sh_get_coeff((pol_sol+os+j),l,m,a_or_b,convert_to_dt,value);
- fprintf(out,"%11.4e ",value[0]);
+ sh_get_coeff((pol_sol+os+j),l,m,
+ a_or_b,convert_to_dt,value);
+ fprintf(out,"%11.4e ",(double)value[0]);
} /* end u_1 .. u_4 nu_1 nu_2 loop */
fprintf(out,"\n");
} /* end layer loop */
@@ -361,7 +376,7 @@
/*
print toroidal solution vector (kernel), not expansion
*/
-void hc_print_toroidal_solution(double *tvec,int lmax,
+void hc_print_toroidal_solution(HC_PREC *tvec,int lmax,
struct hcs *hc,int l_max_out,
char *filename,
hc_boolean verbose)
@@ -382,7 +397,8 @@
for(l=1;l <= ll;l++){
for(os=i=0;i < nl;i++,os+=lmaxp1)
fprintf(out,"%3i %16.7e %16.7e %16.7e\n",
- l,hc->r[i],tvec[os+l],tvec[os2+os+l]);
+ l,(double)hc->r[i],(double)tvec[os+l],
+ (double)tvec[os2+os+l]);
fprintf(out,"\n");
}
fclose(out);
@@ -395,10 +411,10 @@
*/
-void hc_print_vtk(FILE *out,float *xloc,float *xvec,
+void hc_print_vtk(FILE *out,HC_PREC *xloc,HC_PREC *xvec,
int npoints_orig,int nlay,
hc_boolean binary,int shps_d,
- float *xscalar,int nlon, int nlat)
+ HC_PREC *xscalar,int nlon, int nlat)
{
int i,ilay,ndata,poff,j,ndata_d,nele_lay,nele_x,nele_y,
npe,npe1,ncon[12],k,nleft,nlon_m1,npoints,
@@ -406,7 +422,7 @@
nele_brick = 8,nele_tri = 6;
hc_boolean little_endian;
- float xtmp[3],r,spole[3],npole[3];
+ HC_PREC xtmp[3],r,spole[3],npole[3];
/* determine machine type */
little_endian = hc_is_little_endian();
/* */
@@ -434,7 +450,7 @@
hc_print_be_float((xloc+poff),3,out,little_endian);
else
fprintf(out,"%.6e %.6e %.6e\n",
- xloc[poff],xloc[poff+1],xloc[poff+2]);
+ (double)xloc[poff],(double)xloc[poff+1],(double)xloc[poff+2]);
}
/*
south and north poles, add two, to go to npoints per layer
@@ -448,13 +464,13 @@
if(binary)
hc_print_be_float((xtmp),3,out,little_endian);
else
- fprintf(out,"%.6e %.6e %.6e\n",xtmp[0],xtmp[1],xtmp[2]);
+ fprintf(out,"%.6e %.6e %.6e\n",(double)xtmp[0],(double)xtmp[1],(double)xtmp[2]);
/* north pole */
xtmp[2] = r;
if(binary)
hc_print_be_float((xtmp),3,out,little_endian);
else
- fprintf(out,"%.6e %.6e %.6e\n",xtmp[0],xtmp[1],xtmp[2]);
+ fprintf(out,"%.6e %.6e %.6e\n",(double)xtmp[0],(double)xtmp[1],(double)xtmp[2]);
}
/* */
nele_x = nlon;nlon_m1=nlon-1; /* elements in longitude */
@@ -577,7 +593,7 @@
if(binary)
hc_print_be_float((xscalar+poff),1,out,little_endian);
else{
- fprintf(out,"%.6e ",xscalar[poff]);
+ fprintf(out,"%.6e ",(double)xscalar[poff]);
if(i%20 == 0)fprintf(out,"\n");
}
if(i < nlon)
@@ -585,11 +601,11 @@
if((i >= tl) && (i < tr))
npole[0] += xscalar[poff];
}
- spole[0] /= (float)nlon;
- npole[0] /= (float)nlon;
+ spole[0] /= (HC_PREC)nlon;
+ npole[0] /= (HC_PREC)nlon;
if(!binary){ /* ascii */
fprintf(out,"\n");
- fprintf(out,"%.6e %.6e\n",spole[0],npole[0]);
+ fprintf(out,"%.6e %.6e\n",(double)spole[0],(double)npole[0]);
}else{ /* binary */
hc_print_be_float(spole,1,out,little_endian);
hc_print_be_float(npole,1,out,little_endian);
@@ -614,12 +630,12 @@
if(binary) /* binary */
hc_print_be_float((xvec+poff),3,out,little_endian);
else /* ascii */
- fprintf(out,"%.6e %.6e %.6e\n",xvec[poff],xvec[poff+1],
- xvec[poff+2]);
+ fprintf(out,"%.6e %.6e %.6e\n",(double)xvec[poff],(double)xvec[poff+1],
+ (double)xvec[poff+2]);
}
for(k=0;k<3;k++){
- spole[k] /= (float)nlon;
- npole[k] /= (float)nlon;
+ spole[k] /= (HC_PREC)nlon;
+ npole[k] /= (HC_PREC)nlon;
}
if(binary){
/* binary */
@@ -627,40 +643,72 @@
hc_print_be_float(npole,3,out,little_endian);
}else{
/* ascii */
- fprintf(out,"%.6e %.6e %.6e\n",spole[0],spole[1],spole[2]);
- fprintf(out,"%.6e %.6e %.6e\n",npole[0],npole[1],npole[2]);
+ fprintf(out,"%.6e %.6e %.6e\n",(double)spole[0],(double)spole[1],(double)spole[2]);
+ fprintf(out,"%.6e %.6e %.6e\n",(double)npole[0],(double)npole[1],(double)npole[2]);
}
}
}
/*
print big endian binary to file, no matter what hardware
-
+ in HC_BIN_PREC precision
*/
-void hc_print_be_float(float *x, int n, FILE *out,
+int hc_print_be_float(HC_PREC *x, int n, FILE *out,
hc_boolean little_endian)
{
- int i;
- float *xcopy;
- static size_t len = sizeof(float);
+ int i,ret;
+ HC_BIN_PREC *xcopy;
+ const size_t len = sizeof(HC_BIN_PREC);
+ hc_svecalloc(&xcopy,n,"hc_print_be_float");
+ for(i=0;i<n;i++) /* have to make copy */
+ xcopy[i] = (HC_BIN_PREC)x[i];
+
if(little_endian){
/* need to flip the byte order */
- hc_svecalloc(&xcopy,n,"hc_print_be_float");
- memcpy(xcopy,x,len*n);
for(i=0;i < n;i++)
hc_flip_byte_order((void *)(xcopy+i),len);
- fwrite(xcopy,len,n,out);
- free(xcopy);
+ ret= fwrite(xcopy,len,n,out);
}else{
/* can write as is */
- fwrite(x,len,n,out);
+ ret = fwrite(xcopy,len,n,out);
}
+ free(xcopy);
+ return ret;
}
+
+/* print binary to file */
+int hc_print_float(HC_PREC *x, int n, FILE *out)
+{
+ int i,ret;
+ HC_BIN_PREC *xcopy;
+ const size_t len = sizeof(HC_BIN_PREC);
+ hc_svecalloc(&xcopy,n,"hc_print_float");
+ for(i=0;i<n;i++) /* have to make copy */
+ xcopy[i] = (HC_BIN_PREC)x[i];
+ ret = fwrite(xcopy,len,n,out);
+ free(xcopy);
+ return ret;
+}
+/* read binary from file */
+int hc_read_float(HC_PREC *x, int n, FILE *in)
+{
+ int i,ret;
+ HC_BIN_PREC *xcopy;
+ const size_t len = sizeof(HC_BIN_PREC);
+ hc_svecalloc(&xcopy,n,"hc_read_float");
+ ret = fread(xcopy,len,n,in);
+ for(i=0;i<ret;i++) /* have to make copy */
+ x[i] = (HC_PREC)xcopy[i];
+ free(xcopy);
+ return ret;
+}
+
+
void hc_print_be_int(int *x, int n, FILE *out,
hc_boolean little_endian)
{
int i, *xcopy;
- static size_t len = sizeof(int);
+ const size_t len = sizeof(int);
if(little_endian){
/* need to flip the byte order */
hc_ivecalloc(&xcopy,n,"hc_print_be_int");
@@ -736,9 +784,11 @@
sh_print_parameters_to_file((exp+2),1,i,hc->nradp2,HC_Z_DEPTH(hc->r[i]),out,FALSE,binary,verbose);
sh_print_coefficients_to_file((exp+2),1,out,fac,binary,verbose);
if(verbose>2)fprintf(stderr,"hc_print_dens_anom: z: %8.3f (f1: %6.3f f2: %6.3f) %3i/%3i pow: %10.3e %10.3e %10.3e\n",
- HC_Z_DEPTH(hc->r[i]),f1,f2,i+1,hc->nradp2,
- sqrt(sh_total_power((hc->dens_anom+i1))),sqrt(sh_total_power((hc->dens_anom+i2))),
- sqrt(sh_total_power((exp+2))));
+ (double)HC_Z_DEPTH(hc->r[i]),
+ (double)f1,(double)f2,i+1,hc->nradp2,
+ (double)sqrt(sh_total_power((hc->dens_anom+i1))),
+ (double)sqrt(sh_total_power((hc->dens_anom+i2))),
+ (double)sqrt(sh_total_power((exp+2))));
}
sh_free_expansion(exp,3);
Modified: mc/1D/hc/trunk/hc_polsol.c
===================================================================
--- mc/1D/hc/trunk/hc_polsol.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/hc_polsol.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -194,9 +194,10 @@
prop_s1,prop_s2,nvisp1,nzero,n6,ninho,nl=0,ip1;
int newprp,newpot,jpb,inho2,ibv,indx[3],a_or_b,ilayer,lmax,
nprops_max,jsol;
- HC_PREC *xprem;
- double *b,du1,du2,el,rnext,drho,dadd;
- double amat[3][3],bvec[3],u[4],poten[2],unew[4],potnew[2],clm[2];
+ double *xprem;
+ HC_HIGH_PREC *b,du1,du2,el,rnext,drho,dadd;
+ HC_HIGH_PREC amat[3][3],bvec[3],u[4],poten[2],
+ unew[4],potnew[2],clm[2];
/*
structures which hold u[6][4] type arrays
*/
@@ -260,15 +261,15 @@
*/
- hc_dvecalloc(&hc->props,prop_s1 * lmax,"hc_polsol");
- hc_dvecalloc(&hc->ppots,prop_s2 * lmax,"hc_polsol");
+ hc_hvecalloc(&hc->props,prop_s1 * lmax,"hc_polsol");
+ hc_hvecalloc(&hc->ppots,prop_s2 * lmax,"hc_polsol");
}
}else{
/*
propagator recomputed and reallocated each time
*/
- hc_dvecalloc(&hc->props,prop_s1,"hc_polsol");
- hc_dvecalloc(&hc->ppots,prop_s2,"hc_polsol");
+ hc_hvecalloc(&hc->props,prop_s1,"hc_polsol");
+ hc_hvecalloc(&hc->ppots,prop_s2,"hc_polsol");
}
if(!hc->psp.abg_init){
//
@@ -280,7 +281,7 @@
hc->psp.beta = -4.0 * HC_PI * (hc->g*1e3) * (hc->re*1e2) / hc->gacc;
if(verbose)
fprintf(stderr,"hc_polsol: alpha: %.8f beta: %.8f\n",
- hc->psp.alpha,hc->psp.beta);
+ (double)hc->psp.alpha,(double)hc->psp.beta);
/*
geoid scaling factor hc->gacc shoud be grav[nprops] for
@@ -315,9 +316,9 @@
arrays that go with nprops
*/
- hc_dvecalloc(&hc->rprops,nprops_max,"hc_polsol: rprop");
- hc_dvecalloc(&hc->pvisc,nprops_max,"hc_polsol");
- hc_dvecalloc(&hc->den,nprops_max,"hc_polsol");
+ hc_hvecalloc(&hc->rprops,nprops_max,"hc_polsol: rprop");
+ hc_hvecalloc(&hc->pvisc,nprops_max,"hc_polsol");
+ hc_hvecalloc(&hc->den,nprops_max,"hc_polsol");
/* initialize qwrite with zeroes! */
hc->qwrite = (hc_boolean *)calloc(nprops_max,sizeof(hc_boolean));
if(!hc->qwrite)
@@ -433,8 +434,8 @@
if(fabs(hc->den[i]) > HC_EPS_PREC)
i2++;
fprintf(stderr,"hc_polsol: prop: i: %3i(%3i) r: %8.5f v: %8.3f den: %12g ninho: %3i/%3i\n",
- i+1,hc->nprops,hc->rprops[i],
- hc->pvisc[i],hc->den[i],i2,inho);
+ i+1,hc->nprops,(double)hc->rprops[i],
+ (double)hc->pvisc[i],(double)hc->den[i],i2,inho);
}
}
if(!hc->psp.rho_init){
@@ -456,10 +457,11 @@
*/
if(!hc->prem_init)
HC_ERROR("hc_polsol","PREM wasn't initialized for compressible");
- hc_vecalloc(&xprem,hc->prem->np,"hc_polsol: rho");
+ hc_dvecalloc(&xprem,hc->prem->np,"hc_polsol: rho");
for(i=0;i < hc->nprops+1;i++){
- ilayer = prem_find_layer_x(hc->rprops[i],1.0,
- hc->prem->r,10,hc->prem->np,
+ ilayer = prem_find_layer_x((double)hc->rprops[i],1.0,
+ hc->prem->r,
+ 10,hc->prem->np,
xprem);
hc->rho_zero[i] = prem_compute_pval(xprem,
(hc->prem->crho+ilayer*hc->prem->np),
@@ -496,7 +498,7 @@
if(verbose >= 3)
for(i=0;i < hc->nprops+2;i++)
fprintf(stderr,"i: %3i nprops: %3i r(i): %11g rho: %11g\n",
- i,hc->nprops,hc->rprops[i],hc->rho_zero[i]);
+ i,hc->nprops,(double)hc->rprops[i],(double)hc->rho_zero[i]);
//
// begin l loop
Modified: mc/1D/hc/trunk/hc_propagator.c
===================================================================
--- mc/1D/hc/trunk/hc_propagator.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/hc_propagator.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -7,7 +7,7 @@
*/
#include "hc.h"
-void hc_evalpa(int l,double r1,double r2,double visc, double *p)
+void hc_evalpa(int l,HC_HIGH_PREC r1,HC_HIGH_PREC r2,HC_HIGH_PREC visc, HC_HIGH_PREC *p)
{
//
// ****************************************************************
@@ -18,7 +18,7 @@
// * INTERFACE WITH THE EXISTING PROGRAMS. *
// ****************************************************************
//
- double den1,den2,f[4],r,rlm1,rlp1,rmlm2,rml,v2,rs;
+ HC_HIGH_PREC den1,den2,f[4],r,rlm1,rlp1,rmlm2,rml,v2,rs;
long int np[4][4][4];
int lp1,lp2,lp3,lm1,lm2,lpp,lmm,l2p3,l2p1,l2m1,lltp1,lltp2;
int i,j,k,os1,os2;
@@ -32,8 +32,8 @@
// OTHER VARIABLES:
// NP: INTEGER FACTORS (FUNCTIONS OF L) IN THE EXPRESSIONS
// FOR THE ELEMENTS OF P,
- // F: DOUBLE PRECISION FACTORS (FUNCTIONS OF R AND L) FOR P,
- // R IS THE RADIUS RATIO R/R0 (DOUBLE PRECISION).
+ // F: HC_HIGH_PREC PRECISION FACTORS (FUNCTIONS OF R AND L) FOR P,
+ // R IS THE RADIUS RATIO R/R0 (HC_HIGH_PREC PRECISION).
//
r = r2 / r1;
@@ -56,13 +56,13 @@
v2 = visc * 2.0;
rs = r * r;
- rlm1 = pow(r,(double)lm1);
+ rlm1 = pow(r,(HC_HIGH_PREC)lm1);
rlp1 = rlm1 * rs;
- rmlm2 = pow(r,(double)(-lp2));
+ rmlm2 = pow(r,(HC_HIGH_PREC)(-lp2));
rml = rmlm2 * rs;
- den1 = (double)(l2p1 * l2p3);
- den2 = (double)(l2p1 * l2m1);
+ den1 = (HC_HIGH_PREC)(l2p1 * l2p3);
+ den2 = (HC_HIGH_PREC)(l2p1 * l2m1);
f[0] = rlp1 / den1;
f[1] = rlm1 / den2;
@@ -120,7 +120,7 @@
os2 = os1 + j;
p[os2] = 0.0;
for(k=0;k < 4;k++)
- p[os2] += ((double)(np[i][j][k])) * f[k];
+ p[os2] += ((HC_HIGH_PREC)(np[i][j][k])) * f[k];
}
}
@@ -137,7 +137,7 @@
p[3*4+1] *= v2;
}
-void hc_evppot(int l,double ratio, double *ppot)
+void hc_evppot(int l,HC_HIGH_PREC ratio, HC_HIGH_PREC *ppot)
{
// ********************************************
// * THIS SUBROUTINE OBTAINS THE POTENTIALS *
@@ -145,19 +145,19 @@
// * RATIO = R1/R2. *
// ********************************************
//
- double c,expf1,expf2,x,xp1;
+ HC_HIGH_PREC c,expf1,expf2,x,xp1;
// PASSED PARAMETERS: L: DEGREE,
// RATIO: R1 / R2 (R2>R1),
// PPOT: THE POTENTIALS PROPAGATOR [4]
//
- // DOUBLE PRECISION: C: COEFFICIENT MULTIPLYING PROPAGATOR,
+ // HC_HIGH_PREC PRECISION: C: COEFFICIENT MULTIPLYING PROPAGATOR,
// EXPF1: RATIO**L, EXPONENTIAL FACTOR IN PPOT,
// EXPF2: (1/RATIO)**(L+1), EXP. FACTOR IN PPOT,
// X: DEGREE (L),
// XP1: X+1.
- x = (double)l;
+ x = (HC_HIGH_PREC)l;
c = 1.0 / (2.0 * x + 1.0);
xp1 = x + 1.0;
expf1 = pow(ratio,x);
Modified: mc/1D/hc/trunk/hc_solve.c
===================================================================
--- mc/1D/hc/trunk/hc_solve.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/hc_solve.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -54,7 +54,7 @@
(if set to FALSE, can
compare with Benhard's
densub densol output */
- double *tvec;
+ HC_PREC *tvec;
if(!hc->initialized)
HC_ERROR("hc_solve","hc structure not initialized");
@@ -287,7 +287,7 @@
data has to be initialized, eg. as NULL
*/
void hc_compute_sol_spatial(struct hcs *hc, struct sh_lms *sol_w,
- float **sol_x, hc_boolean verbose)
+ HC_PREC **sol_x, hc_boolean verbose)
{
int i,i3,np,np2,np3,os;
static int ntype = 3;
@@ -298,7 +298,7 @@
expansions for r, pol, tor
*/
/* allocate space for spatial solution*/
- hc_svecrealloc(sol_x,np3*hc->nradp2,"sol_x");
+ hc_vecrealloc(sol_x,np3*hc->nradp2,"sol_x");
/*
compute the plm factors
*/
Modified: mc/1D/hc/trunk/hc_torsol.c
===================================================================
--- mc/1D/hc/trunk/hc_torsol.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/hc_torsol.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -51,7 +51,7 @@
void hc_torsol(struct hcs *hc,
int nrad,int nvis,int lmax,HC_PREC *r,
HC_PREC **rv,HC_PREC **visc, struct sh_lms *pvel_tor,
- struct sh_lms *tor_sol,double *tvec,
+ struct sh_lms *tor_sol,HC_HIGH_PREC *tvec,
hc_boolean verbose)
{
//
@@ -61,8 +61,8 @@
// * FIRST ELEMENT AT THE SURFACE IS 1.0. *
// ****************************************************************
//
- double coef,*vecnor,hold,rlast,rnext,tloc[2],*tvec1,*tvec2;
- double exp_fac[2],p[2][2],diflog,el,elp2,elm1,efdiff;
+ HC_HIGH_PREC coef,*vecnor,hold,rlast,rnext,tloc[2],*tvec1,*tvec2;
+ HC_HIGH_PREC exp_fac[2],p[2][2],diflog,el,elp2,elm1,efdiff;
int l,jvisp1,jvis,i,j,nvisp1,lmaxp1,os;
hc_boolean qvis;
//
@@ -228,7 +228,7 @@
// set tvec(l,nradp2-1,0) = 1.0 and normalize all vectors to
// this
//
- hc_dvecalloc(&vecnor,lmaxp1,"hc_torsol: vecnor");
+ hc_hvecalloc(&vecnor,lmaxp1,"hc_torsol: vecnor");
os = (hc->nradp2-1) * lmaxp1;
vecnor[0] = 1.0;
for(l=1;l < lmaxp1;l++)
Modified: mc/1D/hc/trunk/main.c
===================================================================
--- mc/1D/hc/trunk/main.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/main.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -57,7 +57,7 @@
FILE *out;
struct hc_parameters p[1]; /* parameters */
char filename[HC_CHAR_LENGTH],file_prefix[10];
- float *sol_spatial = NULL; /* spatial solution,
+ HC_PREC *sol_spatial = NULL; /* spatial solution,
e.g. velocities */
HC_PREC corr[2]; /* correlations */
HC_PREC vl[4][3],v[4],dv; /* for viscosity scans */
@@ -227,7 +227,7 @@
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]);
+ fprintf(stdout,"%10.7f %10.7f\n",(double)corr[0],(double)corr[1]);
}else{
/*
print geoid solution
Modified: mc/1D/hc/trunk/prem_util.c
===================================================================
--- mc/1D/hc/trunk/prem_util.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/prem_util.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -147,7 +147,8 @@
/*
get density in SI units [kg/m^3] at non-dim radius rnd
*/
-void prem_get_rho(double *rho,double rnd, struct prem_model *prem)
+void prem_get_rho(double *rho,double rnd,
+ struct prem_model *prem)
{
int ilayer,os;
double *x;
Modified: mc/1D/hc/trunk/rick_fft_c.c
===================================================================
--- mc/1D/hc/trunk/rick_fft_c.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/rick_fft_c.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -68,11 +68,11 @@
void rick_realft_nr(SH_RICK_PREC *rdata, int n, int isign)
{
SH_RICK_PREC c1,c2,h1r,h1i,h2r,h2i;
- double theta,wi,wpi,wpr,wr,wtemp;
+ SH_RICK_HIGH_PREC theta,wi,wpi,wpr,wr,wtemp;
int i,n2p3,ilim,i1,i2,i3,i4,n2;
static int negunity = -1,unity = 1;
- theta = RICK_PI/(double)(n);
+ theta = RICK_PI/(SH_RICK_HIGH_PREC)(n);
wr = 1.0;
wi = 0.0;
@@ -145,9 +145,9 @@
// +-NN/2 to -1 in the standard fashion.
//
// local
- double tempr,tempi;
- // this should be double precision locally, regardless
- double wr,wi,wpr,wpi,wtemp,theta,temp2;
+ SH_RICK_HIGH_PREC tempr,tempi;
+ // this should be SH_RICK_HIGH_PREC precision locally, regardless
+ SH_RICK_HIGH_PREC wr,wi,wpr,wpi,wtemp,theta,temp2;
int n,m,i,j,mmax,istep;
Modified: mc/1D/hc/trunk/rick_sh_c.c
===================================================================
--- mc/1D/hc/trunk/rick_sh_c.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/rick_sh_c.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -6,8 +6,8 @@
// ivec is set to 1
// */
-void rick_compute_allplm(int lmax,int ivec,float *plm,
- float *dplm, struct rick_module *rick)
+void rick_compute_allplm(int lmax,int ivec,SH_RICK_PREC *plm,
+ SH_RICK_PREC *dplm, struct rick_module *rick)
{
int i,os;
@@ -29,7 +29,7 @@
void rick_compute_allplm_reg(int lmax,int ivec,SH_RICK_PREC *plm,
SH_RICK_PREC *dplm, struct rick_module *rick,
- float *theta, int ntheta)
+ SH_RICK_PREC *theta, int ntheta)
{
int i,os;
@@ -108,9 +108,9 @@
exit(-1);
}
/* allocate memory */
- hc_svecalloc(&plm,rick->nlat*rick->lmsize,"rick_shc2d: mem 1");
+ rick_vecalloc(&plm,rick->nlat*rick->lmsize,"rick_shc2d: mem 1");
if(ivec)
- hc_svecalloc(&dplm,rick->nlat*rick->lmsize,"rick_shc2d: mem 2");
+ rick_vecalloc(&dplm,rick->nlat*rick->lmsize,"rick_shc2d: mem 2");
//
// compute the Plm first
rick_compute_allplm(lmax,ivec,plm,dplm,rick);
@@ -132,8 +132,8 @@
void rick_shc2d_reg(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,
+ struct rick_module *rick,SH_RICK_PREC *theta, int ntheta,
+ SH_RICK_PREC *phi,int nphi,
hc_boolean save_sincos_fac)
{
SH_RICK_PREC *plm,*dplm;
@@ -142,9 +142,9 @@
exit(-1);
}
/* allocate memory */
- hc_svecalloc(&plm,ntheta*rick->lmsize,"rick_shc2d_reg: mem 1");
+ rick_vecalloc(&plm,ntheta*rick->lmsize,"rick_shc2d_reg: mem 1");
if(ivec)
- hc_svecalloc(&dplm,ntheta*rick->lmsize,"rick_shc2d_reg: mem 2");
+ rick_vecalloc(&dplm,ntheta*rick->lmsize,"rick_shc2d_reg: mem 2");
//
// compute the Plm first for all theta values
rick_compute_allplm_reg(lmax,ivec,plm,dplm,rick,theta,ntheta);
@@ -166,7 +166,7 @@
// */
void rick_shc2d_pre(SH_RICK_PREC *cslm,SH_RICK_PREC *dslm,
int lmax,SH_RICK_PREC *plm, SH_RICK_PREC *dplm,
- int ivec,float *rdatax,float *rdatay,
+ int ivec,SH_RICK_PREC *rdatax,SH_RICK_PREC *rdatay,
struct rick_module *rick)
{
/* //
@@ -317,16 +317,16 @@
*/
void rick_shc2d_pre_reg(SH_RICK_PREC *cslm,SH_RICK_PREC *dslm,
int lmax,SH_RICK_PREC *plm, SH_RICK_PREC *dplm,
- int ivec,float *rdatax,float *rdatay,
- struct rick_module *rick, float *theta,
- int ntheta,float *phi,int nphi,
+ int ivec,SH_RICK_PREC *rdatax,SH_RICK_PREC *rdatay,
+ struct rick_module *rick, SH_RICK_PREC *theta,
+ int ntheta,SH_RICK_PREC *phi,int nphi,
my_boolean save_sincos_fac)
{
/* //
// Legendre functions are precomputed
// */
- double dpdt,dpdp,mphi,sin_theta;
- float *loc_plma=NULL,*loc_plmb=NULL;
+ SH_RICK_HIGH_PREC dpdt,dpdp,mphi,sin_theta;
+ SH_RICK_PREC *loc_plma=NULL,*loc_plmb=NULL;
int i,j,k,k2,m,ios1,ios2,ios3,l,lmaxp1,lm1,idata;
if(!rick->initialized){
fprintf(stderr,"rick_shc2d_pre_reg: error: initialize modules first\n");
@@ -349,13 +349,13 @@
*/
if((!save_sincos_fac)||(!rick->sin_cos_saved)){
- hc_svecrealloc(&rick->sfac,nphi*lmaxp1,"rick_shc2d_pre_reg");
- hc_svecrealloc(&rick->cfac,nphi*lmaxp1,"rick_shc2d_pre_reg");
+ rick_vecrealloc(&rick->sfac,nphi*lmaxp1,"rick_shc2d_pre_reg");
+ rick_vecrealloc(&rick->cfac,nphi*lmaxp1,"rick_shc2d_pre_reg");
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);
+ mphi = (SH_RICK_HIGH_PREC)m*(SH_RICK_HIGH_PREC)phi[i];
+ rick->sfac[ios1] = (SH_RICK_PREC)sin(mphi);
+ rick->cfac[ios1] = (SH_RICK_PREC)cos(mphi);
}
}
if(save_sincos_fac)
@@ -367,8 +367,8 @@
scalar
*/
- hc_svecrealloc(&loc_plma,ntheta*rick->lmsize,"rick_shc2d_pre_reg 3");
- hc_svecrealloc(&loc_plmb,ntheta*rick->lmsize,"rick_shc2d_pre_reg 4");
+ rick_vecrealloc(&loc_plma,ntheta*rick->lmsize,"rick_shc2d_pre_reg 3");
+ rick_vecrealloc(&loc_plmb,ntheta*rick->lmsize,"rick_shc2d_pre_reg 4");
for(i=ios1=0;i < ntheta;i++){ /* theta dependent array */
for(k=k2=0;k < rick->lmsize;k++,k2+=2,ios1++){
loc_plma[ios1] = cslm[k2 ] * plm[ios1];
@@ -450,21 +450,21 @@
}
if(!save_sincos_fac){
- hc_svecrealloc(&rick->sfac,1,"");
- hc_svecrealloc(&rick->cfac,1,"");
+ rick_vecrealloc(&rick->sfac,1,"");
+ rick_vecrealloc(&rick->cfac,1,"");
}
}
/* completely irregular output */
void rick_shc2d_irreg(SH_RICK_PREC *cslm,SH_RICK_PREC *dslm,
- int lmax,int ivec,float *rdatax,float *rdatay,
- struct rick_module *rick, float *theta,
- float *phi,int npoints)
+ int lmax,int ivec,SH_RICK_PREC *rdatax,SH_RICK_PREC *rdatay,
+ struct rick_module *rick, SH_RICK_PREC *theta,
+ SH_RICK_PREC *phi,int npoints)
{
/* //
// Legendre functions are precomputed
// */
- double dpdt,dpdp,mphi,sin_theta,sfac,cfac;
- float *plm=NULL,*dplm=NULL;
+ SH_RICK_HIGH_PREC dpdt,dpdp,mphi,sin_theta,sfac,cfac;
+ SH_RICK_PREC *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");
@@ -481,7 +481,7 @@
scalar
*/
- hc_svecrealloc(&plm,rick->lmsize,"rick_shc2d_irreg: mem 1");
+ rick_vecrealloc(&plm,rick->lmsize,"rick_shc2d_irreg: mem 1");
for(i=0;i < npoints;i++){
/* get legendre function values */
rick_plmbar1(plm,dplm,ivec,lmax,cos(theta[i]),rick);
@@ -495,7 +495,7 @@
}
//fprintf(stderr,"%5i %5i %11g %11g\n",l,m,cslm[k2],cslm[k2+1]);
if(m != 0){
- mphi = (double)m*(double)phi[i];
+ mphi = (SH_RICK_HIGH_PREC)m*(SH_RICK_HIGH_PREC)phi[i];
rdatax[i] += cslm[k2] * plm[k] * cos(mphi);
rdatax[i] += cslm[k2+1] * plm[k] * sin(mphi);
}else{
@@ -510,8 +510,8 @@
/*
vector harmonics
*/
- hc_svecrealloc(&plm,rick->lmsize,"rick_shc2d_irreg: mem 1");
- hc_svecrealloc(&dplm,rick->lmsize,"rick_shc2d_irreg: mem 2");
+ rick_vecrealloc(&plm,rick->lmsize,"rick_shc2d_irreg: mem 1");
+ rick_vecrealloc(&dplm,rick->lmsize,"rick_shc2d_irreg: mem 2");
for(i=0;i < npoints;i++){
/* get legendre function values */
rick_plmbar1(plm,dplm,ivec,lmax,cos(theta[i]),rick);
@@ -536,9 +536,9 @@
dpdp *= (SH_RICK_PREC)rick->ell_factor[lm1]; /* d_phi (P_lm) factor */
if(m){
- mphi = (double)m*(double)phi[i];
- sfac = (float)sin(mphi);
- cfac = (float)cos(mphi);
+ mphi = (SH_RICK_HIGH_PREC)m*(SH_RICK_HIGH_PREC)phi[i];
+ sfac = (SH_RICK_PREC)sin(mphi);
+ cfac = (SH_RICK_PREC)cos(mphi);
}else{
mphi = sfac = 0.0;
cfac = 1.0;
@@ -620,8 +620,8 @@
// local
SH_RICK_PREC *plm,*dplm;
/* allocate memory */
- hc_svecalloc(&plm,rick->nlat*rick->lmsize,"rick_shd2c: mem 1");
- hc_svecalloc(&dplm,rick->nlat*rick->lmsize,"rick_shd2c: mem 2");
+ rick_vecalloc(&plm,rick->nlat*rick->lmsize,"rick_shd2c: mem 1");
+ rick_vecalloc(&dplm,rick->nlat*rick->lmsize,"rick_shd2c: mem 2");
// check
if(!rick->initialized){
fprintf(stderr,"rick_shd2c: error: initialize first\n");
@@ -882,11 +882,11 @@
//
// those will be used by plmbar to store some of the factors
//
- hc_svecalloc(&rick->plm_f1,rick->lmsize,"rick_init 4");
- hc_svecalloc(&rick->plm_f2,rick->lmsize,"rick_init 5");
- hc_svecalloc(&rick->plm_fac1,rick->lmsize,"rick_init 6");
- hc_svecalloc(&rick->plm_fac2,rick->lmsize,"rick_init 7");
- hc_svecalloc(&rick->plm_srt,rick->nlon,"rick_init 8");
+ rick_vecalloc(&rick->plm_f1,rick->lmsize,"rick_init 4");
+ rick_vecalloc(&rick->plm_f2,rick->lmsize,"rick_init 5");
+ rick_vecalloc(&rick->plm_fac1,rick->lmsize,"rick_init 6");
+ rick_vecalloc(&rick->plm_fac2,rick->lmsize,"rick_init 7");
+ rick_vecalloc(&rick->plm_srt,rick->nlon,"rick_init 8");
rick->sin_cos_saved = FALSE;
@@ -1013,7 +1013,7 @@
// SH_RICK_PREC,intent(inout), dimension (lmsize) p, dp
//
// local
- double plm,pm1,pm2,pmm,sintsq,fnum,fden;
+ SH_RICK_HIGH_PREC plm,pm1,pm2,pmm,sintsq,fnum,fden;
//
int i,l,m,k,kstart,l2,mstop,lmaxm1;
if(!rick->initialized){
@@ -1131,7 +1131,7 @@
for(l = 2;l<=lmax;l++){ // now all (l,0)
k += l;
l2 = l * 2;
- plm = ((double)(l2-1) * z * pm1 - (double)(l-1) * pm2)/((double)l);
+ plm = ((SH_RICK_HIGH_PREC)(l2-1) * z * pm1 - (SH_RICK_HIGH_PREC)(l-1) * pm2)/((SH_RICK_HIGH_PREC)l);
p[k] = rick->plm_srt[l2] * plm;
pm2 = pm1;
pm1 = plm;
@@ -1212,7 +1212,7 @@
// local variables
//
int i,j,m;
- double p1,p2,p3,pp,xl,xm,z,z1;
+ SH_RICK_HIGH_PREC p1,p2,p3,pp,xl,xm,z,z1;
//
pp = 0.0; /* for copiler */
m=(n+1)/2;
@@ -1227,9 +1227,9 @@
for(j = 1;j <= n;j++){
p3 = p2;
p2 = p1;
- p1 = ((2.0*(double)j - 1.0) * z * p2 -(j- 1.0)*p3)/(double)j;
+ p1 = ((2.0*(SH_RICK_HIGH_PREC)j - 1.0) * z * p2 -(j- 1.0)*p3)/(SH_RICK_HIGH_PREC)j;
}
- pp = (double)n*(z*p1-p2)/(z*z-1.0);
+ pp = (SH_RICK_HIGH_PREC)n*(z*p1-p2)/(z*z-1.0);
z1 = z;
z = z1-p1/pp;
}
Modified: mc/1D/hc/trunk/sh.h
===================================================================
--- mc/1D/hc/trunk/sh.h 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/sh.h 2012-07-16 12:18:58 UTC (rev 20526)
@@ -12,6 +12,8 @@
same for Rick's spherical harmonic routines
*/
+#define SH_RICK_PRECISION HC_PRECISION
+
#include "sh_rick.h"
/*
@@ -116,7 +118,7 @@
*/
my_boolean save_plm;
/* layer indicators if nset != 1 */
- float *z;
+ HC_PREC *z;
/* scalar only? ivec=0 or velocities? ivec=1 */
int ivec;
/* number of expansions per set */
@@ -132,7 +134,7 @@
*/
int tnpoints; /* number of the total datapoints */
/* data structure */
- float *data;
+ HC_PREC *data;
my_boolean spatial_init, initialized;
};
/*
Modified: mc/1D/hc/trunk/sh_ana.c
===================================================================
--- mc/1D/hc/trunk/sh_ana.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/sh_ana.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -55,7 +55,7 @@
struct sh_lms *exp;
hc_boolean verbose = TRUE, use_3d = FALSE, short_format = FALSE,read_grd =FALSE,
binary = FALSE, print_spatial_base = FALSE;
- float *data, zlabel = 0,*flt_dummy;
+ HC_PREC *data, zlabel = 0,*flt_dummy;
struct ggrd_gt ggrd[2];
SH_RICK_PREC *dummy;
HC_PREC fac[3] = {1.,1.,1.};
@@ -150,7 +150,7 @@
/* intialize expansion first */
sh_allocate_and_init(&exp,shps*nset,lmax,type,ivec,verbose,FALSE);
/* make room for data */
- hc_svecalloc(&data,shps * exp->npoints,"sh_ana");
+ hc_vecalloc(&data,shps * exp->npoints,"sh_ana");
if(print_spatial_base){
/*
print out spatial basis
Modified: mc/1D/hc/trunk/sh_corr.c
===================================================================
--- mc/1D/hc/trunk/sh_corr.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/sh_corr.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -125,11 +125,11 @@
if(verbose)
fprintf(stderr,"%s: computing linear correlation per degree from %i to %i \n",argv[0],lmin,llim);
for(l=lmin;l<=llim;l++)
- fprintf(stdout,"%5i %14.7e\n",l,sh_correlation_per_degree(exp1,exp2,l,l));
+ fprintf(stdout,"%5i %14.7e\n",l,(double)sh_correlation_per_degree(exp1,exp2,l,l));
}else{
if(verbose)
fprintf(stderr,"%s: computing total linear correlation from %i to %i\n",argv[0],lmin,llim);
- fprintf(stdout,"%14.7e\n",sh_correlation_per_degree(exp1,exp2,lmin,llim));
+ fprintf(stdout,"%14.7e\n",(double)sh_correlation_per_degree(exp1,exp2,lmin,llim));
}
Modified: mc/1D/hc/trunk/sh_exp.c
===================================================================
--- mc/1D/hc/trunk/sh_exp.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/sh_exp.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -119,11 +119,8 @@
/*
use single precision vector
*/
-#ifdef SH_RICK_DOUBLE_PRECISION
- hc_dvecalloc(&exp->alm,exp->n_lm,"sh_init_expansion");
-#else
- hc_svecalloc(&exp->alm,exp->n_lm,"sh_init_expansion");
-#endif
+ rick_vecalloc(&exp->alm,exp->n_lm,"sh_init_expansion");
+
sh_clear_alm(exp);
/*
@@ -227,10 +224,10 @@
*/
HC_CPREC sh_total_power(struct sh_lms *exp)
{
- float *power;
+ HC_PREC *power;
double sum;
int l;
- hc_svecalloc(&power,exp->lmaxp1,"sh_total_power");
+ hc_vecalloc(&power,exp->lmaxp1,"sh_total_power");
sh_compute_power_per_degree(exp,power);
for(sum=0.0,l=0;l<=exp->lmax;l++)
sum += (double)(2.0*(HC_CPREC)l+1.0) * (double)power[l];
@@ -245,7 +242,7 @@
*/
void sh_compute_power_per_degree(struct sh_lms *exp,
- float *power)
+ HC_PREC *power)
{
int l,m;
HC_CPREC value[2];
@@ -273,7 +270,8 @@
HC_PREC sh_correlation_per_degree(struct sh_lms *exp1, struct sh_lms *exp2, int lmin,int lmax)
{
int l,m;
- HC_CPREC sum[3],tmp,atmp,btmp,ctmp,dtmp,value1[2],value2[2];
+ HC_CPREC sum[3],tmp,atmp,btmp,ctmp,value1[2],value2[2];
+ double dtmp;
hc_boolean need_b;
sum[0]=sum[1]=sum[2]=0.0;
@@ -337,7 +335,7 @@
FILE *out, hc_boolean short_format,
hc_boolean binary,hc_boolean verbose)
{
- HC_BIN_PREC fz;
+ HC_PREC fz;
/*
print
@@ -345,11 +343,11 @@
*/
if(binary){
- fz = (HC_BIN_PREC)zlabel;
+ fz = (HC_PREC)zlabel;
fwrite(&exp[0].lmax,sizeof(int),1,out);
if(!short_format){
fwrite(&ilayer,sizeof(int),1,out);
- fwrite(&fz,sizeof(HC_BIN_PREC),1,out);
+ hc_print_float(&fz,1,out);
fwrite(&nset,sizeof(int),1,out);
fwrite(&shps,sizeof(int),1,out);
fwrite(&exp[0].type,sizeof(int),1,out);
@@ -357,7 +355,8 @@
}else{
if(!short_format)
fprintf(out,"%6i %6i %.8e %6i %2i %2i ",
- exp[0].lmax,ilayer,zlabel,nset,
+ exp[0].lmax,ilayer,
+ (double)zlabel,nset,
shps,exp[0].type);
else
fprintf(out,"%6i ",
@@ -424,7 +423,7 @@
hc_boolean verbose)
{
int input1[2],input2[3];
- HC_BIN_PREC fz;
+ HC_PREC fz;
double dtmp;
/*
read
@@ -439,7 +438,7 @@
*lmax = input1[0];
}else{
if(fread(input1,sizeof(int),2,in)+
- fread(&fz, sizeof(HC_BIN_PREC),1,in) +
+ hc_read_float(&fz,1,in) +
fread(input2,sizeof(int),3,in) != 6)
return FALSE;
*lmax = input1[0];
@@ -459,7 +458,7 @@
lmax,ilayer,&dtmp,nset,shps,type)!=6){
return FALSE;
}
- *zlabel = (HC_CPREC)dtmp;
+ *zlabel = (HC_PREC)dtmp;
}
}
if(short_format){
@@ -544,8 +543,8 @@
hc_boolean verbose)
{
int j,l,m;
- HC_CPREC value[2];
- HC_BIN_PREC fvalue[2];
+ HC_PREC value[2];
+ HC_PREC fvalue[2];
/*
test other expansions this set
*/
@@ -570,9 +569,9 @@
we are using internally
*/
sh_get_coeff((exp+j),l,m,2,TRUE,value);
- fvalue[0] = (HC_BIN_PREC)value[0]*fac[j];
- fvalue[1] = (HC_BIN_PREC)value[1]*fac[j];
- fwrite(fvalue, sizeof(HC_BIN_PREC), 2, out);
+ fvalue[0] = value[0]*fac[j];
+ fvalue[1] = value[1]*fac[j];
+ hc_print_float(fvalue, 2, out);
}
}else{
for(l=0;l <= exp[0].lmax;l++){
@@ -582,7 +581,8 @@
convention */
sh_get_coeff((exp+j),l,m,2,TRUE,value);
fprintf(out,"%15.7e %15.7e\t",
- value[0]*fac[j],value[1]*fac[j]);
+ (double)(value[0]*fac[j]),
+ (double)(value[1]*fac[j]));
}
fprintf(out,"\n");
} /* end m loop */
@@ -611,7 +611,7 @@
{
int j,k,l,m,lmax_loc;
HC_CPREC value[2]={0,0};
- HC_BIN_PREC fvalue[2]={0,0};
+ HC_PREC fvalue[2]={0,0};
if(lmax < 0)
lmax_loc = exp[0].lmax;
else
@@ -635,7 +635,7 @@
for(l=0;l <= lmax_loc;l++)
for(m=0;m <= l;m++)
for(j=0;j < shps;j++){
- if(fread(fvalue, sizeof(HC_BIN_PREC),2,in)!=2){
+ if(hc_read_float(fvalue,2,in)!=2){
fprintf(stderr,"sh_read_coefficients_from_file: read error: set %i l %i m %i\n",
j+1,l,m);
exit(-1);
@@ -650,9 +650,9 @@
for(l=0;l <= lmax_loc;l++)
for(m=0;m <= l;m++)
for(j=0;j < shps;j++){
- if(fscanf(in,"%lf %lf",value,(value+1))!=2){
+ if(fscanf(in,HC_TWO_FLT_FORMAT,value,(value+1))!=2){
fprintf(stderr,"sh_read_coefficients_from_file: read error: set %i l %i m %i, last val: %g %g\n",
- j+1,l,m,value[0],value[1]);
+ j+1,l,m,(double)value[0],(double)value[1]);
exit(-1);
}
/* read in real, Dahlen & Tromp normalized coefficients and
@@ -688,7 +688,8 @@
for(m=0;m <= l;m++){
sh_get_coeff(exp,l,m,2,FALSE,value);
if(fabs(value[0])+fabs(value[1]) > 1e-8)
- fprintf(out,"%5i %5i %15.7e %15.7e\n",l,m,value[0],value[1]);
+ fprintf(out,"%5i %5i %15.7e %15.7e\n",l,m,
+ (double)value[0],(double)value[1]);
}
}
}
@@ -714,7 +715,8 @@
*/
void sh_read_spatial_data_from_file(struct sh_lms *exp, FILE *in,
my_boolean use_3d,
- int shps, float *data, float *z)
+ int shps, HC_PREC *data,
+ HC_PREC *z)
{
struct ggrd_gt *gdummy=(struct ggrd_gt *)NULL;
sh_read_spatial_data(exp,in,gdummy,FALSE,use_3d,shps,data,z);
@@ -722,7 +724,8 @@
/* similar for grd files */
void sh_read_spatial_data_from_grd(struct sh_lms *exp, struct ggrd_gt *ggrd,
my_boolean use_3d,
- int shps, float *data, float *z)
+ int shps, HC_PREC *data,
+ HC_PREC *z)
{
FILE *fdummy=NULL;
sh_read_spatial_data(exp,fdummy,ggrd,TRUE,use_3d,shps,data,z);
@@ -733,10 +736,11 @@
*/
void sh_read_spatial_data(struct sh_lms *exp, FILE *in, struct ggrd_gt *ggrd,my_boolean use_grd,
- my_boolean use_3d, int shps, float *data, float *z)
+ my_boolean use_3d, int shps,
+ HC_PREC *data, HC_PREC *z)
{
- float lon,lat,xp[3];
+ HC_PREC lon,lat,xp[3];
int j,k;
double dvalue;
/*
@@ -772,7 +776,8 @@
rick_pix2ang(j,exp->lmax,(xp+HC_THETA),(xp+HC_PHI),
&exp->rick);
#else
- rick_f90_pix2ang(&j,&exp->lmax,(xp+HC_THETA),(xp+HC_PHI));
+ rick_f90_pix2ang(&j,&exp->lmax,(SH_RICK_PREC)(xp+HC_THETA),
+ (SH_RICK_PREC)(xp+HC_PHI));
#endif
break;
#ifdef HC_USE_SPHEREPACK
@@ -792,7 +797,7 @@
/*
read in lon lat
*/
- if(fscanf(in,"%f %f",&lon,&lat) != 2){
+ if(fscanf(in,HC_TWO_FLT_FORMAT,&lon,&lat) != 2){
fprintf(stderr,"sh_read_spatial_data: error: lon lat format: pixel %i: read error\n",
(int)j);
exit(-1);
@@ -803,7 +808,7 @@
/*
read in lon lat z[i]
*/
- if(fscanf(in,"%f %f %f",&lon,&lat,z) != 3){
+ if(fscanf(in,HC_THREE_FLT_FORMAT,&lon,&lat,z) != 3){
fprintf(stderr,"sh_read_spatial_data: error: lon lat z format: pixel %i: read error\n",
(int)j);
exit(-1);
@@ -820,7 +825,7 @@
for(k=0;k < shps;k++){
if(!ggrd_grdtrack_interpolate_tp((double)xp[HC_THETA],(double)xp[HC_PHI],(ggrd+k),&dvalue,FALSE,FALSE)){
fprintf(stderr,"sh_read_spatial_data: interpolation error grd %i, lon %g lat %g\n",
- k+1,PHI2LON(xp[HC_PHI]),THETA2LAT(xp[HC_THETA]));
+ k+1,(double)PHI2LON(xp[HC_PHI]),(double)THETA2LAT(xp[HC_THETA]));
exit(-1);
}
data[k*exp[0].npoints+j] = dvalue;
@@ -830,7 +835,7 @@
read data
*/
for(k=0;k < shps;k++){
- if(fscanf(in,"%f",(data+k*exp[0].npoints+j))!=1){
+ if(fscanf(in,HC_FLT_FORMAT,(data+k*exp[0].npoints+j))!=1){
fprintf(stderr,"sh_read_spatial_data: error: scalar format: pixel %i: read error\n",
(int)j);
exit(-1);
@@ -849,7 +854,7 @@
fprintf(stderr,"sh_read_model_spatial_data: error: pixel %i coordinate mismatch:\n",
(int)j);
fprintf(stderr,"sh_read_model_spatial_data: orig: %g, %g file: %g, %g\n",
- PHI2LON(xp[HC_PHI]),THETA2LAT(xp[HC_THETA]),lon,lat);
+ (double)PHI2LON(xp[HC_PHI]),(double)THETA2LAT(xp[HC_THETA]),lon,lat);
exit(-1);
}
}
@@ -880,14 +885,14 @@
*/
void sh_compute_spatial_basis(struct sh_lms *exp, FILE *out,
hc_boolean use_3d,
- float z, float **x,
+ HC_PREC z, HC_PREC **x,
int out_mode, /* 0: file 1: store */
hc_boolean verbose)
{
int j,os,inc;
- float xp[3];
+ HC_PREC xp[3];
if(out_mode) /* make room for storing x,y,z */
- hc_svecrealloc(x,exp->npoints*(2+((use_3d)?(1):(0))),"sh_compute_spatial_basis");
+ hc_vecrealloc(x,exp->npoints*(2+((use_3d)?(1):(0))),"sh_compute_spatial_basis");
inc = (use_3d)?(3):(2);
for(j=os=0;j < exp->npoints;j++,os+=inc){
/*
@@ -935,23 +940,23 @@
if(out_mode){ /* save */
if(!use_3d){
/* save in lon lat format */
- *(*x+os) = (float) PHI2LON(xp[HC_PHI]);
- *(*x+os+1) = (float) THETA2LAT(xp[HC_THETA]);
+ *(*x+os) = (HC_PREC) PHI2LON(xp[HC_PHI]);
+ *(*x+os+1) = (HC_PREC) THETA2LAT(xp[HC_THETA]);
}else{
/* save in lon lat z format */
- *(*x+os) = (float) PHI2LON(xp[HC_PHI]);
- *(*x+os+1) = (float) THETA2LAT(xp[HC_THETA]);
- *(*x+os+2) = (float) z;
+ *(*x+os) = (HC_PREC) PHI2LON(xp[HC_PHI]);
+ *(*x+os+1) = (HC_PREC) THETA2LAT(xp[HC_THETA]);
+ *(*x+os+2) = (HC_PREC) z;
}
}else{ /* write to file */
if(!use_3d){
/* write in lon lat format */
- fprintf(out,"%14.7f %14.7f\n",PHI2LON(xp[HC_PHI]),
- THETA2LAT(xp[HC_THETA]));
+ fprintf(out,"%14.7f %14.7f\n",(double)PHI2LON(xp[HC_PHI]),
+ (double)THETA2LAT(xp[HC_THETA]));
}else{
/* write in lon lat z format */
- fprintf(out,"%14.7f %14.7f %g\n",PHI2LON(xp[HC_PHI]),
- THETA2LAT(xp[HC_THETA]),z);
+ fprintf(out,"%14.7f %14.7f %g\n",(double)PHI2LON(xp[HC_PHI]),
+ (double)THETA2LAT(xp[HC_THETA]),(double)z);
}
}
}
@@ -990,8 +995,8 @@
IVEC = 1
*/
-void sh_compute_spectral(float *data, int ivec,
- hc_boolean save_plm,float **plm,
+void sh_compute_spectral(HC_PREC *data, int ivec,
+ hc_boolean save_plm,HC_PREC **plm,
struct sh_lms *exp, hc_boolean verbose)
{
if(save_plm){
@@ -1089,7 +1094,7 @@
*/
void sh_compute_spatial(struct sh_lms *exp, int ivec,
hc_boolean save_plm,SH_RICK_PREC **plm,
- float *data, hc_boolean verbose)
+ HC_PREC *data, hc_boolean verbose)
{
if((!exp[0].spectral_init)||(ivec && !exp[1].spectral_init)){
fprintf(stderr,"sh_compute_spatial: coefficients set not initialized, ivec: %i\n",
@@ -1166,10 +1171,11 @@
*/
void sh_compute_spatial_reg(struct sh_lms *exp, int ivec,
- hc_boolean save_plm,SH_RICK_PREC **plm,
- float *theta, int ntheta, float *phi,int nphi,
- float *data, hc_boolean verbose,
- hc_boolean save_sincos_fac)
+ hc_boolean save_plm,SH_RICK_PREC **plm,
+ HC_PREC *theta, int ntheta,
+ HC_PREC *phi,int nphi,
+ HC_PREC *data, hc_boolean verbose,
+ hc_boolean save_sincos_fac)
{
int npoints;
npoints = nphi * ntheta;
@@ -1221,8 +1227,8 @@
}
}
void sh_compute_spatial_irreg(struct sh_lms *exp, int ivec,
- float *theta, float *phi,int npoints,
- float *data, hc_boolean verbose)
+ HC_PREC *theta, HC_PREC *phi,int npoints,
+ HC_PREC *data, hc_boolean verbose)
{
if((!exp[0].spectral_init)||(ivec && !exp[1].spectral_init)){
fprintf(stderr,"sh_compute_spatial_irreg: coefficients set not initialized, ivec: %i\n",
@@ -1282,7 +1288,7 @@
jlim=(ivec)?(3):(1);
for(i=0;i < n_plm;i++){ /* number of points loop */
for(j=0;j < jlim;j++) /* scalar or scalar + pol? */
- fprintf(out,"%16.7e ",plm[j*n_plm+i]);
+ fprintf(out,"%16.7e ",(double)plm[j*n_plm+i]);
fprintf(out,"\n");
}
break;
@@ -1291,7 +1297,7 @@
jlim=(ivec)?(2):(1);
for(i=0;i < n_plm;i++){ /* number of points loop */
for(j=0;j < jlim;j++) /* scalar or scalar + pol? */
- fprintf(out,"%16.7e ",plm[j*n_plm+i]);
+ fprintf(out,"%16.7e ",(double)plm[j*n_plm+i]);
fprintf(out,"\n");
}
break;
@@ -1321,11 +1327,11 @@
*/
void sh_print_spatial_data_to_file(struct sh_lms *exp, int shps,
- float *data, hc_boolean use_3d,
- float z, FILE *out)
+ HC_PREC *data, hc_boolean use_3d,
+ HC_PREC z, FILE *out)
{
int j,k;
- float lon,lat;
+ HC_PREC lon,lat;
for(j=0;j < exp[0].npoints;j++){
/*
get coordinates
@@ -1334,21 +1340,21 @@
/* print coordinates */
if(!use_3d){
/* print lon lat */
- fprintf(out,"%14.7f %14.7f\t",lon,lat);
+ fprintf(out,"%14.7f %14.7f\t",(double)lon,(double)lat);
}else{
/* print lon lat z[i] */
- fprintf(out,"%14.7f %14.7f %14.7f\t",lon,lat,z);
+ fprintf(out,"%14.7f %14.7f %14.7f\t",(double)lon,(double)lat,(double)z);
}
for(k=0;k < shps;k++) /* loop through all scalars */
- fprintf(out,"%14.7e ",data[j+exp[0].npoints*k]);
+ fprintf(out,"%14.7e ",(double)data[j+exp[0].npoints*k]);
fprintf(out,"\n");
} /* end points in lateral space loop */
}
void sh_get_coordinates(struct sh_lms *exp,
- int i, float *lon, float *lat)
+ int i, HC_PREC *lon, HC_PREC *lat)
{
- float xp[3];
+ HC_PREC xp[3];
switch(exp->type){
#ifdef HC_USE_HEALPIX
@@ -1398,13 +1404,13 @@
*/
void sh_print_reg_spatial_data_to_file(struct sh_lms *exp, int shps,
- float *data, hc_boolean use_3d,
- float z, float *theta,int ntheta,
- float *phi,int nphi,
+ HC_PREC *data, hc_boolean use_3d,
+ HC_PREC z, HC_PREC *theta,int ntheta,
+ HC_PREC *phi,int nphi,
FILE *out)
{
int i,j,k,l,npoints;
- float lon,lat;
+ HC_PREC lon,lat;
npoints = nphi * ntheta;
/*
get coordinates
@@ -1424,13 +1430,13 @@
/* print coordinates */
if(!use_3d){
/* print lon lat */
- fprintf(out,"%12.3f %12.3f\t",lon,lat);
+ fprintf(out,"%12.3f %12.3f\t",(double)lon,(double)lat);
}else{
/* print lon lat z[i] */
- fprintf(out,"%12.3f %12.3f %14.7f\t",lon,lat,z);
+ fprintf(out,"%12.3f %12.3f %14.7f\t",(double)lon,(double)lat,(double)z);
}
for(k=0;k<shps;k++) /* loop through all scalars */
- fprintf(out,"%14.7e ",data[l+npoints*k]);
+ fprintf(out,"%14.7e ",(double)data[l+npoints*k]);
fprintf(out,"\n");
}
}
@@ -1453,12 +1459,12 @@
*/
void sh_print_irreg_spatial_data_to_file(struct sh_lms *exp, int shps,
- float *data, hc_boolean use_3d,
- float z, float *theta,float *phi,int npoints,
+ HC_PREC *data, hc_boolean use_3d,
+ HC_PREC z, HC_PREC *theta,HC_PREC *phi,int npoints,
FILE *out)
{
int i,k;
- float lon,lat;
+ HC_PREC lon,lat;
/*
get coordinates
*/
@@ -1475,13 +1481,13 @@
/* print coordinates */
if(!use_3d){
/* print lon lat */
- fprintf(out,"%14.7f %14.7f\t",lon,lat);
+ fprintf(out,"%14.7f %14.7f\t",(double)lon,(double)lat);
}else{
/* print lon lat z[i] */
- fprintf(out,"%14.7f %14.7f %14.7f\t",lon,lat,z);
+ fprintf(out,"%14.7f %14.7f %14.7f\t",(double)lon,(double)lat,(double)z);
}
for(k=0;k < shps;k++) /* loop through all scalars */
- fprintf(out,"%14.7e ",data[i+npoints*k]);
+ fprintf(out,"%14.7e ",(double)data[i+npoints*k]);
fprintf(out,"\n");
}
break;
@@ -1510,7 +1516,8 @@
plm: will be re-allocated, has to be passed at least as NULL
*/
-void sh_compute_plm(struct sh_lms *exp,int ivec,SH_RICK_PREC **plm,
+void sh_compute_plm(struct sh_lms *exp,int ivec,
+ SH_RICK_PREC **plm,
hc_boolean verbose)
{
if(!exp->plm_computed){
@@ -1523,7 +1530,7 @@
/*
allocate
*/
- hc_svecrealloc(plm,exp->tn_plm,"sh_compute_plm");
+ rick_vecalloc(plm,exp->tn_plm,"sh_compute_plm");
/*
compute the Legendre polynomials
*/
@@ -1602,8 +1609,10 @@
plm: will be re-allocated, has to be passed at least as NULL
*/
-void sh_compute_plm_reg(struct sh_lms *exp,int ivec,SH_RICK_PREC **plm,
- hc_boolean verbose, float *theta, int npoints)
+void sh_compute_plm_reg(struct sh_lms *exp,int ivec,
+ SH_RICK_PREC **plm,
+ hc_boolean verbose,
+ HC_PREC *theta, int npoints)
{
/* */
exp->tn_plm_irr = (1+ivec) * exp->lmsmall2 * npoints;
@@ -1619,7 +1628,7 @@
allocate
*/
exp->old_tnplm_irr = exp->tn_plm_irr;
- hc_svecrealloc(plm,exp->old_tnplm_irr,"sh_compute_plm_reg");
+ rick_vecrealloc(plm,exp->old_tnplm_irr,"sh_compute_plm_reg");
/*
compute the Legendre polynomials
*/
Modified: mc/1D/hc/trunk/sh_model.c
===================================================================
--- mc/1D/hc/trunk/sh_model.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/sh_model.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -96,7 +96,7 @@
/*
layer indicator attributes
*/
- model->z = (float *)calloc(model->nset,sizeof(float));
+ model->z = (HC_PREC *)calloc(model->nset,sizeof(HC_PREC));
if(!model->z)
HC_MEMERROR("sh_init_model: z");
/*
@@ -196,7 +196,7 @@
hc_boolean verbose)
{
int i;
- float **flt_dummy=NULL;
+ HC_PREC **flt_dummy=NULL;
if(verbose)
fprintf(stderr,"sh_print_model_spatial_basis: printing spatial basis for nset %i expansions\n",
model->nset);
@@ -229,7 +229,7 @@
*/
void sh_read_model_spatial_data(struct sh_lms_model *model,
- float **data,FILE *in,
+ HC_PREC **data,FILE *in,
hc_boolean verbose)
{
int i,j;
@@ -242,8 +242,8 @@
/*
allocate space for all the data points
*/
- hc_svecrealloc(data, model->tnpoints,
- "sh_read_model_spatial_data");
+ hc_vecrealloc(data, model->tnpoints,
+ "sh_read_model_spatial_data");
/*
check all sets, number of points have to be the same for each
expansion
@@ -261,7 +261,8 @@
}
for(i=0;i < model->nset;i++)
/* read in the data for this layer */
- sh_read_spatial_data_from_file((model->exp+i*model->shps),in,(model->nset!=1)?(TRUE):(FALSE),
+ sh_read_spatial_data_from_file((model->exp+i*model->shps),
+ in,(model->nset!=1)?(TRUE):(FALSE),
model->shps,
(*data+model->shps*model->exp[i*model->shps].npoints),
(model->z+i));
@@ -274,7 +275,7 @@
*/
void sh_compute_model_spectral(struct sh_lms_model *model,
- float *data,
+ HC_PREC *data,
hc_boolean verbose)
{
int i;
@@ -316,7 +317,7 @@
*/
void sh_compute_model_spatial(struct sh_lms_model *model,
- float **data,hc_boolean verbose)
+ HC_PREC **data,hc_boolean verbose)
{
int i;
const int unity = 1, zero = 0;
@@ -329,7 +330,7 @@
/*
make room for data
*/
- hc_svecrealloc(data,model->tnpoints,"sh_compute_model_spatial");
+ hc_vecrealloc(data,model->tnpoints,"sh_compute_model_spatial");
if(((model->ivec)&&(model->shps != 3))||
((!model->ivec)&&(model->shps != 1))){
fprintf(stderr,"sh_compute_model_spatial: error: ivec: %i nshp: %i\n",
@@ -364,7 +365,7 @@
*/
void sh_print_model_spatial_data(struct sh_lms_model *model,
- float *data, FILE *out,
+ HC_PREC *data, FILE *out,
hc_boolean verbose)
{
int i,j,os;
Modified: mc/1D/hc/trunk/sh_power.c
===================================================================
--- mc/1D/hc/trunk/sh_power.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/sh_power.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -36,7 +36,7 @@
{
int type,lmax,shps,ilayer,nset,ivec,i,l;
hc_boolean verbose = TRUE, short_format = FALSE ,binary = FALSE;
- float *power = NULL;
+ HC_PREC *power = NULL;
HC_PREC fac[3] = {1.,1.,1.},zlabel;
struct sh_lms *exp;
if(argc>1){
@@ -76,7 +76,7 @@
sh_allocate_and_init(&exp,shps,lmax, type,ivec,verbose,FALSE);
sh_read_coefficients_from_file(exp,shps,-1,stdin,binary,fac,verbose);
/* get space */
- hc_svecrealloc(&power,exp->lmaxp1 * shps,"sh_power");
+ hc_vecrealloc(&power,exp->lmaxp1 * shps,"sh_power");
/* compute the powers */
for(i=0;i<shps;i++)
sh_compute_power_per_degree((exp+i),(power+i*exp->lmaxp1));
Modified: mc/1D/hc/trunk/sh_rick.h
===================================================================
--- mc/1D/hc/trunk/sh_rick.h 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/sh_rick.h 2012-07-16 12:18:58 UTC (rev 20526)
@@ -32,15 +32,31 @@
//
-#ifdef SH_RICK_DOUBLE_PRECISION
+#if SH_RICK_PRECISION == 32
+
+#define SH_RICK_HIGH_PREC long double
+#define SH_RICK_PREC long double
+#define rick_vecalloc hc_vecalloc
+#define rick_vecrealloc hc_vecrealloc
+#define SH_RICK_FLT_FMT "%Lf"
+
+#elif SH_RICK_PRECISION == 16 /* double */
+
+#define SH_RICK_HIGH_PREC double
#define SH_RICK_PREC double
#define rick_vecalloc hc_dvecalloc
+#define rick_vecrealloc hc_dvecrealloc
#define SH_RICK_FLT_FMT "%lf"
-#else
+
+#else /* single */
+
+#define SH_RICK_HIGH_PREC double
#define SH_RICK_PREC float
#define rick_vecalloc hc_svecalloc
+#define rick_vecrealloc hc_svecrealloc
#define SH_RICK_FLT_FMT "%f"
+
#endif
#ifndef my_boolean
Modified: mc/1D/hc/trunk/sh_syn.c
===================================================================
--- mc/1D/hc/trunk/sh_syn.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/sh_syn.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -23,12 +23,12 @@
*/
hc_boolean verbose = TRUE, short_format = FALSE ,short_format_ivec = FALSE ,binary = FALSE;
int regular_basis = 0;
- double w,e,s,n,dx,dy;
+ HC_PREC w,e,s,n,dx,dy;
/* */
FILE *in;
- float *data,*theta,*phi;
+ HC_PREC *data,*theta,*phi;
/* spacing for reg_ular output */
- double dphi,x,y,dtheta;
+ HC_PREC dphi,x,y,dtheta;
HC_PREC fac[3] = {1.,1.,1.},zlabel;
SH_RICK_PREC *dummy;
struct sh_lms *exp;
@@ -49,27 +49,27 @@
short_format_ivec = TRUE;
}
if(argc > 3){
- sscanf(argv[3],"%lf",&w);
+ sscanf(argv[3],HC_FLT_FORMAT,&w);
if(w == 999)
regular_basis = -1;
else
regular_basis = 1;
}
if(argc > 4)
- sscanf(argv[4],"%lf",&e);
+ sscanf(argv[4],HC_FLT_FORMAT,&e);
if(argc > 5)
- sscanf(argv[5],"%lf",&s);
+ sscanf(argv[5],HC_FLT_FORMAT,&s);
if(argc > 6)
- sscanf(argv[6],"%lf",&n);
+ sscanf(argv[6],HC_FLT_FORMAT,&n);
if(argc > 7)
- sscanf(argv[7],"%lf",&dx);
+ sscanf(argv[7],HC_FLT_FORMAT,&dx);
if(argc > 8)
- sscanf(argv[8],"%lf",&dy);
+ sscanf(argv[8],HC_FLT_FORMAT,&dy);
else
dy = dx;
if((argc > 9)|| (argc < 0)){
fprintf(stderr,"usage: %s [short_format, %i] [short_ivec, %i] [w, %g] [e, %g] [s, %g] [n, %g] [dx, %g] [dy, dx] (in that order)\n",
- argv[0],short_format,short_format_ivec,w,e,s,n,dx);
+ argv[0],short_format,short_format_ivec,(double)w,(double)e,(double)s,(double)n,(double)dx);
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");
@@ -105,7 +105,7 @@
*/
if(verbose)
fprintf(stderr,"sh_syn: using regular spaced grid with -R%g/%g/%g/%g -I%g/%g spacing\n",
- w,e,s,n,dx,dy);
+ (double)w,(double)e,(double)s,(double)n,(double)dx,(double)dy);
if((w > e)||(s>n)||(s<-90)||(s>90)||(n<-90)||(n>90)){
fprintf(stderr,"%s: range error\n",argv[0]);
exit(-1);
@@ -115,7 +115,7 @@
s += dy/2;
n -= dy/2;
fprintf(stderr,"sh_syn: vector fields: adjusting to -R%g/%g/%g/%g\n",
- w,e,s,n);
+ (double)w,(double)e,(double)s,(double)n);
}
/* */
dphi = DEG2RAD(dx);
@@ -125,13 +125,13 @@
npoints = nphi * ntheta;
/* */
- hc_svecalloc(&phi,nphi,"sh_shsyn");
- hc_svecalloc(&theta,ntheta,"sh_shsyn");
+ hc_vecalloc(&phi,nphi,"sh_shsyn");
+ hc_vecalloc(&theta,ntheta,"sh_shsyn");
for(x=LON2PHI(w),i=0;i < nphi;i++,x += dphi)
phi[i] = x;
for(y = LAT2THETA(n),j=0;j < ntheta;y += dtheta,j++)
theta[j] = y;
- hc_svecalloc(&data,npoints * shps,"sh_shsyn data");
+ hc_vecalloc(&data,npoints * shps,"sh_shsyn data");
/* compute the expansion */
sh_compute_spatial_reg(exp,ivec,FALSE,&dummy,
theta,ntheta,phi,nphi,data,
@@ -146,30 +146,30 @@
if(verbose)
fprintf(stderr,"sh_syn: reading locations lon lat from stdin for expansion\n");
npoints = 0;
- hc_svecalloc(&phi,1,"sh_syn");
- hc_svecalloc(&theta,1,"sh_syn");
+ hc_vecalloc(&phi,1,"sh_syn");
+ hc_vecalloc(&theta,1,"sh_syn");
in = fopen("tmp.lonlat","r");
if(!in){
fprintf(stderr,"sh_syn: error, could not open tmp.lonlat for reading lon lat locations\n");
exit(-1);
}
- while(fscanf(in,"%lf %lf",&dphi,&dtheta)==2){
+ while(fscanf(in,HC_TWO_FLT_FORMAT,&dphi,&dtheta)==2){
phi[npoints] = LON2PHI(dphi);
theta[npoints] = LAT2THETA(dtheta);
npoints++;
- hc_svecrealloc(&phi,npoints+1,"sh_syn");
- hc_svecrealloc(&theta,npoints+1,"sh_syn");
+ hc_vecrealloc(&phi,npoints+1,"sh_syn");
+ hc_vecrealloc(&theta,npoints+1,"sh_syn");
}
if(verbose)
fprintf(stderr,"sh_syn: read %i locations lon lat from tmp.lonlat for expansion\n",npoints);
fclose(in);
- hc_svecalloc(&data,npoints * shps,"sh_shsyn data");
+ hc_vecalloc(&data,npoints * shps,"sh_shsyn data");
sh_compute_spatial_irreg(exp,ivec,theta,phi,npoints,data,verbose);
sh_print_irreg_spatial_data_to_file(exp,shps,data,(nset>1)?(TRUE):(FALSE),
zlabel,theta,phi,npoints,stdout);
}else{ /* use the built in spatial basis (Gaussian) */
/* expansion */
- hc_svecalloc(&data,exp[0].npoints * shps,"sh_syn");
+ hc_vecalloc(&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),
Modified: mc/1D/hc/trunk/simple_test.c
===================================================================
--- mc/1D/hc/trunk/simple_test.c 2012-07-15 15:31:06 UTC (rev 20525)
+++ mc/1D/hc/trunk/simple_test.c 2012-07-16 12:18:58 UTC (rev 20526)
@@ -2,7 +2,6 @@
#include <math.h>
-void hc_dvecalloc(double **,int, char *);
#define vshd2c vshd2c_
#define gauleg gauleg_
@@ -54,10 +53,4 @@
0.,0.,-cpol[i*2],-cpol[i*2+1],-ctor[i*2],-ctor[i*2+1]);
}
-void hc_dvecalloc(double **x,int n,char *message)
-{
- *x = (double *)malloc(sizeof(double)*(size_t)n);
- if(! (*x)){
- fprintf(stderr,"%s: memory error\n",message);
- }
-}
+
More information about the CIG-COMMITS
mailing list