[cig-commits] r18779 - in mc/1D/hc/trunk: . anisotropic_viscosity
becker at geodynamics.org
becker at geodynamics.org
Sat Jul 16 15:37:33 PDT 2011
Author: becker
Date: 2011-07-16 15:37:33 -0700 (Sat, 16 Jul 2011)
New Revision: 18779
Modified:
mc/1D/hc/trunk/Makefile
mc/1D/hc/trunk/anisotropic_viscosity/Anisotropic_viscosity.c
mc/1D/hc/trunk/ggrd_grdtrack_util.c
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/hc_extract_spatial.c
mc/1D/hc/trunk/hc_misc.c
mc/1D/hc/trunk/hc_output.c
mc/1D/hc/trunk/sh_exp.c
mc/1D/hc/trunk/sh_syn.c
Log:
Added basic VTK output functionality and hc_extract_spatial
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2011-07-16 19:19:23 UTC (rev 18778)
+++ mc/1D/hc/trunk/Makefile 2011-07-16 22:37:33 UTC (rev 18779)
@@ -95,7 +95,7 @@
#
HC_SOURCES = sh_exp.c sh_model.c hc_init.c hc_solve.c hc_propagator.c \
hc_polsol.c hc_matrix.c hc_torsol.c hc_output.c hc_input.c \
- hc_misc.c hc_extract_sh_layer.c
+ hc_misc.c hc_extract_sh_layer.c hc_extract_spatial.c
# all C sources
C_SOURCES = $(HC_SOURCES) $(RICK_SRCS) $(GGRD_SRCS)
@@ -140,7 +140,8 @@
sh_tools: $(BDIR)/sh_syn $(BDIR)/sh_corr $(BDIR)/sh_ana $(BDIR)/sh_power \
$(BDIR)/sh_extract_layer
-hc_tools: $(BDIR)/hc $(BDIR)/hc_extract_sh_layer $(BDIR)/rotvec2vel $(BDIR)/print_gauss_lat
+hc_tools: $(BDIR)/hc $(BDIR)/hc_extract_sh_layer $(BDIR)/hc_extract_spatial \
+ $(BDIR)/rotvec2vel $(BDIR)/print_gauss_lat
weird_tools: $(BDIR)/convert_bernhard_dens
@@ -232,6 +233,12 @@
-o $(BDIR)/hc_extract_sh_layer \
-lhc -lrick $(HEAL_LIBS_LINKLINE) $(GGRD_LIBS_LINKLINE) -lm $(LDFLAGS)
+$(BDIR)/hc_extract_spatial: $(LIBS) $(INCS) $(PREM_OBJS) $(ODIR)/hc_extract_spatial.o
+ $(CC) $(LIB_FLAGS) $(ODIR)/hc_extract_spatial.o $(PREM_OBJS) \
+ -o $(BDIR)/hc_extract_spatial \
+ -lhc -lrick $(HEAL_LIBS_LINKLINE) $(GGRD_LIBS_LINKLINE) -lm $(LDFLAGS)
+
+
#
# C function prototyper, strip out GMT version dependent things,
# those are handled in other header
Modified: mc/1D/hc/trunk/anisotropic_viscosity/Anisotropic_viscosity.c
===================================================================
--- mc/1D/hc/trunk/anisotropic_viscosity/Anisotropic_viscosity.c 2011-07-16 19:19:23 UTC (rev 18778)
+++ mc/1D/hc/trunk/anisotropic_viscosity/Anisotropic_viscosity.c 2011-07-16 22:37:33 UTC (rev 18779)
@@ -100,26 +100,28 @@
/*
-transversely isotropic viscosity following Han and Wahr
+transversely isotropic viscosity following Han and Wahr (see their eq. 5)
-\nu_1 = isotropic viscosity, applies for e_31, e_23
-\nu_2 = weak anisotropy, applies for e_31, e_32
-\eta_1 = normal viscosity, (\eta_1+2\nu_1) control e_11, e_22
-\eta_2 = normal viscosity, (\eta_2+2\nu_2) = 2\eta_1 + 2\nu_1, controls e_33
+\nu_1 = strong (isotropic) viscosity, scales e_31, e_11, e_22
+\nu_2 = weak anisotropy, scales for e_31, e_32
+\eta_1 = normal viscosity, (\eta_1+2\nu_1) scales e_11, e_22
+\eta_2 = normal viscosity, (\eta_2+2\nu_2) = 2\eta_1 + 2\nu_1, scales e_33
we use (for consistency with anisotropic viscosity)
-Delta = 1-\nu_2/\nu_1
+\delta = 1-\nu_2/\nu_1
and
\Gamma, such that \eta_1 = \Gamma \nu_1
-\nu_1 is the reference, isotropic viscosity, set to unity here, i.e.
+\nu_1 is the reference, isotropic viscosity, set to unity here,
-\nu_2 = 1 - \Delta ; \eta_1 = \Gamma ; (\eta_2 = 2 (\Gamma-\Delta)); for isotropy \Delta = 0, \Gamma = 0
+for \nu_1 = 1 -->
+\nu_2 = 1 - \delta ; \eta_1 = \Gamma ; \eta_2 = 2\eta_1 + 2\nu_1 - 2\nu_2 = 2\Gamma -2 \delta; isotropy \Delta = 0, \Gamma = 0
+
n[3] is the cartesian direction into which the weak shear points
(ie. routine will rotate the 3 axis into the n direction) and will
normalize n, if not already normalized
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2011-07-16 19:19:23 UTC (rev 18778)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2011-07-16 22:37:33 UTC (rev 18779)
@@ -19,8 +19,8 @@
#endif
#include <math.h>
#include <string.h>
+#include <math.h>
-
void ggrd_init_master(struct ggrd_master *ggrd)
{
ggrd->mat_control = ggrd->mat_control_init = 0;
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2011-07-16 19:19:23 UTC (rev 18778)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2011-07-16 22:37:33 UTC (rev 18779)
@@ -34,6 +34,7 @@
void ggrd_weights(double, double *, int, int, double [(5 +1)][(1 +1)]);
/* grdinttester.c */
/* hc_extract_sh_layer.c */
+/* hc_extract_spatial.c */
/* hc_init.c */
void hc_init_parameters(struct hc_parameters *);
void hc_struc_init(struct hcs **);
@@ -79,6 +80,9 @@
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 *);
/* 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);
@@ -93,6 +97,11 @@
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);
+void hc_print_be_float(float *, int, FILE *, unsigned short);
+unsigned short hc_is_little_endian(void);
+void hc_flip_byte_order(void *, size_t);
+void hc_flipit(void *, void *, size_t);
/* hc_polsol.c */
void hc_polsol(struct hcs *, int, double *, int, double *, unsigned short, struct sh_lms *, unsigned short, int, double *, double *, unsigned short, struct sh_lms *, struct sh_lms *, unsigned short, struct sh_lms *, unsigned short, unsigned short);
/* hc_propagator.c */
@@ -177,6 +186,7 @@
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);
Modified: mc/1D/hc/trunk/hc_extract_spatial.c
===================================================================
--- mc/1D/hc/trunk/hc_extract_spatial.c 2011-07-16 19:19:23 UTC (rev 18778)
+++ mc/1D/hc/trunk/hc_extract_spatial.c 2011-07-16 22:37:33 UTC (rev 18779)
@@ -2,22 +2,24 @@
/*
-extract part of a spherical harmonics solution of a HC run
+extract part of a solution of a HC run and convert to spatial
-$Id$
-
*/
int main(int argc, char **argv)
{
- int ilayer,nsol,i,mode,shps,nset=1,loop,i1,i2;
+ int ilayer,nsol,mode,shps,loop,i1,i2,ivec,lc,ndata,npoints,i,j,
+ poff,ilayer3;
FILE *in;
struct sh_lms *sol=NULL;
struct hcs *model;
- HC_PREC fac[3] = {1.0,1.0,1.0};
- hc_boolean binary = TRUE, verbose = FALSE, short_format = FALSE;
+ HC_PREC zlabel;
+ hc_boolean binary_in = TRUE, verbose = FALSE;
+ float *data,*plm=NULL,*xpos,*xvec,lon,lat,theta,phi,xtmp[3],
+ pvec[3];
+ double polar_base[9];
hc_struc_init(&model);
/*
deal with parameters
@@ -32,40 +34,31 @@
sscanf(argv[2],"%i",&ilayer);
sscanf(argv[3],"%i",&mode);
break;
- case 5:
- sscanf(argv[2],"%i",&ilayer);
- sscanf(argv[3],"%i",&mode);
- sscanf(argv[4],"%i",&i);
- if(i)
- short_format = TRUE;
- else
- short_format = FALSE;
- break;
default:
- fprintf(stderr,"%s: usage\n%s sol.file layer [mode,%i] [short_format, %i]\n\n",
- argv[0],argv[0],mode,short_format);
- fprintf(stderr,"extracts spherical harmonic solution x (vel or str) from HC run\n");
+ fprintf(stderr,"%s: usage\n%s sol.file layer [mode,%i] \n\n",
+ argv[0],argv[0],mode);
+ fprintf(stderr,"extracts spatial solution x (velocity or stress) from HC run\n");
fprintf(stderr,"layer: 1...nset\n");
fprintf(stderr,"\tif ilayer= 1..nset, will print one layer\n");
- fprintf(stderr,"\t -1, will select nset\n");
+ fprintf(stderr,"\t -1, will select nset (the top layer)\n");
fprintf(stderr,"\t -2, will print all layers\n");
- fprintf(stderr,"mode: 1...6\n");
- fprintf(stderr,"\tif mode = 1, will print x_r \n");
- fprintf(stderr,"\t 2, will print x_pol x_tor \n");
- fprintf(stderr,"\t 3, will print x_r x_pol x_tor\n");
+ fprintf(stderr,"mode: 1...4\n");
+ fprintf(stderr,"\tif mode = 1, will print lon lat z v_r \n");
+ fprintf(stderr,"\t 2, will print lon lat z v_theta v_phi \n");
+ fprintf(stderr,"\t 3, will print lon lat z v_r v_theta v_phi\n");
fprintf(stderr,"\t 4, will print the depth levels of all layers\n");
- fprintf(stderr,"\t 5, will print x_pol\n");
- fprintf(stderr,"\t 6, will print x_tor\n");
+ fprintf(stderr,"\t 5, compute all depth levels and write VTK file, ASCII\n");
+ fprintf(stderr,"\t 6, compute all depth levels and write VTK file, BINARY\n");
exit(-1);
break;
}
- if(mode == 4)
+ if((mode == 4)||(mode==5)||(mode==6))
ilayer = -2;
/*
read in solution
*/
in = ggrd_open(argv[1],"r","hc_extract_sh_layer");
- hc_read_sh_solution(model,&sol,in,binary,verbose);
+ hc_read_sh_solution(model,&sol,in,binary_in,verbose);
fclose(in);
nsol = model->nradp2 * 3;
/*
@@ -74,81 +67,75 @@
loop = 0;
if(ilayer == -1)
ilayer = model->nradp2;
- if(ilayer == -2){
+ else if(ilayer == -2){
ilayer = model->nradp2;
loop =1;
}
- if((ilayer<1)||(ilayer > model->nradp2)){
+ if((ilayer < 1)||(ilayer > model->nradp2)){
fprintf(stderr,"%s: ilayer (%i) out of range, use 1 ... %i\n",
argv[0],ilayer,model->nradp2);
exit(-1);
}
+ /* set up layer bounds */
if(loop){
i1=0;i2=model->nradp2-1;
- if(short_format)
- printf("%i\n",model->nradp2);
}else{
i1=ilayer-1;i2 = i1;
}
/* detect number of expansions */
- if((mode == 1)||(mode == 5)||(mode == 6))
+ if(mode == 1){
+ shps = 1; /* r */
+ }else if(mode == 2){
+ shps = 2; /* theta,phi */
+ }else if((mode == 3)||(mode == 5)||(mode==6)){
+ shps = 3; /* r,theta,phi */
+ }else{
shps = 1;
- else if(mode == 2)
- shps = 2;
- else if(mode == 3)
- shps = 3;
- for(ilayer=i1;ilayer <= i2;ilayer++){
+ }
+ /* room for spatial expansion */
+ npoints = (sol+i1*3)->npoints;
+ ndata = npoints * shps;
+ if((mode == 5)||(mode==6)) /* save all layers */
+ hc_svecalloc(&data,model->nradp2 * ndata,"hc_extract_spatial");
+ else
+ hc_svecalloc(&data,ndata,"hc_extract_spatial");
+ for(lc=0,ilayer=i1,ilayer3=ilayer*3;ilayer <= i2;ilayer++,lc++,ilayer3+=3){
/*
output
*/
- if(mode != 4){
- /* SH header */
- if(short_format && loop)
- fprintf(stdout,"%g\n",HC_Z_DEPTH(model->r[ilayer]));
- sh_print_parameters_to_file((sol+ilayer*3),shps,
- ilayer,nset,(HC_PREC)(HC_Z_DEPTH(model->r[ilayer])),
- stdout,short_format,FALSE,verbose);
- }
+
+ zlabel = HC_Z_DEPTH(model->r[ilayer]);
switch(mode){
case 1:
/* */
if(verbose)
- fprintf(stderr,"%s: printing x_r SHE at layer %i (depth: %g)\n",
- argv[0],ilayer,HC_Z_DEPTH(model->r[ilayer]));
- sh_print_coefficients_to_file((sol+ilayer*3),shps,stdout,fac,FALSE,verbose);
+ fprintf(stderr,"%s: printing v_r at layer %i (depth: %g)\n",argv[0],ilayer,zlabel);
+
+ ivec=FALSE;sh_compute_spatial((sol+ilayer3),ivec,TRUE,&plm,data,verbose);
+ sh_print_spatial_data_to_file((sol+ilayer3),shps,data,TRUE,zlabel,stdout);
break;
case 2:
/* */
if(verbose)
- fprintf(stderr,"%s: printing x_pol x_tor SHE at layer %i (depth: %g)\n",
- argv[0],ilayer,HC_Z_DEPTH(model->r[ilayer]));
- sh_print_coefficients_to_file((sol+ilayer*3+1),shps,stdout,fac,FALSE,verbose);
+ fprintf(stderr,"%s: printing v_theta v_phi SHE at layer %i (depth: %g)\n",argv[0],ilayer,zlabel);
+ ivec=TRUE;sh_compute_spatial((sol+ilayer3+1),ivec,TRUE,&plm,data,verbose);
+ sh_print_spatial_data_to_file((sol+ilayer3+1),shps,data,TRUE,zlabel,stdout);
break;
case 3:
- /* mode == 3 */
if(verbose)
- fprintf(stderr,"%s: printing x_r x_pol x_tor SHE at layer %i (depth: %g)\n",
- argv[0],ilayer,HC_Z_DEPTH(model->r[ilayer]));
- sh_print_coefficients_to_file((sol+ilayer*3),shps,stdout,fac,FALSE,verbose);
+ fprintf(stderr,"%s: printing v_r v_theta v_phi SHE at layer %i (depth: %g)\n",argv[0],ilayer,zlabel);
+ ivec=FALSE;sh_compute_spatial((sol+ilayer3), ivec,TRUE,&plm,data,verbose); /* radial */
+ ivec=TRUE; sh_compute_spatial((sol+ilayer3+1),ivec,TRUE,&plm,(data+npoints),verbose); /* theta,phi */
+ sh_print_spatial_data_to_file((sol+ilayer3),shps,data,TRUE,zlabel,stdout);
break;
case 4:
fprintf(stdout,"%5i %11g\n",ilayer,HC_Z_DEPTH(model->r[ilayer]));
break;
- case 5:
- /* */
- if(verbose)
- fprintf(stderr,"%s: printing x_pol SHE at layer %i (depth: %g)\n",
- argv[0],ilayer,HC_Z_DEPTH(model->r[ilayer]));
- sh_print_coefficients_to_file((sol+ilayer*3+1),shps,stdout,fac,FALSE,verbose);
- break;
+ case 5: /* compute all and store */
case 6:
- /* */
- if(verbose)
- fprintf(stderr,"%s: printing x_tor SHE at layer %i (depth: %g)\n",
- argv[0],ilayer,HC_Z_DEPTH(model->r[ilayer]));
- sh_print_coefficients_to_file((sol+ilayer*3+2),shps,stdout,fac,FALSE,verbose);
+ ivec=FALSE;sh_compute_spatial((sol+ilayer3), ivec,TRUE,&plm,(data+lc*ndata),verbose); /* radial */
+ ivec=TRUE; sh_compute_spatial((sol+ilayer3+1),ivec,TRUE,&plm,(data+lc*ndata+npoints),verbose); /* theta,phi */
break;
-
default:
fprintf(stderr,"%s: error, mode %i undefined\n",argv[0],mode);
exit(-1);
@@ -157,6 +144,45 @@
}
/* clear and exit */
sh_free_expansion(sol,nsol);
+ free(plm);
+ /* */
+ if((mode == 5)||(mode==6)){
+ 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"); /* ndata = npoints * 3 */
+ hc_svecalloc(&xvec,model->nradp2 * ndata,"hc_extract_spatial");
+ for(i=0;i < npoints;i++){ /* loop through all points */
+ /* lon lat coordinates */
+ sh_get_coordinates((sol+i1*3),i,&lon,&lat);
+ theta = LAT2THETA(lat);phi = LON2PHI(lon);
+ xtmp[0] = xtmp[1] = sin(theta);
+ xtmp[0] *= cos(phi); /* x */
+ xtmp[1] *= sin(phi); /* y */
+ xtmp[2] = cos(theta); /* z */
+ /* for conversion */
+ calc_polar_base_at_theta_phi(theta,phi,polar_base);
+ for(ilayer=0;ilayer < model->nradp2;ilayer++){
+ /* this is the slow data storage loop but it avoids
+ recomputing the polar basis vector */
+ poff = ilayer * ndata + i*shps; /* point offset */
+ for(j=0;j<3;j++){
+ xpos[poff+j] = xtmp[j] * model->r[ilayer]; /* cartesian coordinates */
+ }
+ /* data are stored a bit weirdly, this makes for lots of
+ jumping around in memory ... */
+ pvec[0] = data[ilayer*ndata+i];
+ pvec[1] = data[ilayer*ndata+npoints+i];
+ pvec[2] = data[ilayer*ndata+npoints*2+i];
+ lonlatpv2cv_with_base(pvec,polar_base,(xvec+poff));
+ }
+ }
+ free(data);
+ /* print in VTK format */
+ hc_print_vtk(stdout,xpos,xvec,npoints,model->nradp2,(mode==6));
+ free(xvec);free(xpos);
+ }else{
+ free(data);
+ }
return 0;
}
Modified: mc/1D/hc/trunk/hc_misc.c
===================================================================
--- mc/1D/hc/trunk/hc_misc.c 2011-07-16 19:19:23 UTC (rev 18778)
+++ mc/1D/hc/trunk/hc_misc.c 2011-07-16 22:37:33 UTC (rev 18779)
@@ -217,3 +217,52 @@
exit(-1);
}
}
+
+/*
+ convert polar vector in r,theta,phi format to cartesian
+ vector x
+
+*/
+
+void lonlatpv2cv(float lon, float lat,
+ float *polar_vec,float *cart_vec)
+{
+ double polar_base[9];
+ float theta,phi;
+ theta = LAT2THETA((double)lat);
+ phi = LON2PHI((double)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)
+{
+ int i;
+ // convert vector
+ for(i=0;i<3;i++){
+ cart_vec[i] = polar_base[i] * polar_vec[0]; /* r,theta,phi */
+ cart_vec[i] += polar_base[3+i] * polar_vec[1];
+ cart_vec[i] += polar_base[6+i] * polar_vec[2];
+ }
+}
+void calc_polar_base_at_theta_phi(float theta, float phi,
+ double *polar_base)
+{
+ double cp,sp,ct,st;
+ // base vecs
+ ct=cos((double)theta);cp=cos((double)phi);
+ st=sin((double)theta);sp=sin((double)phi);
+ //
+ polar_base[0]= st * cp;
+ polar_base[1]= st * sp;
+ polar_base[2]= ct;
+ //
+ polar_base[3]= ct * cp;
+ polar_base[4]= ct * sp;
+ polar_base[5]= -st;
+ //
+ polar_base[6]= -sp;
+ polar_base[7]= cp;
+ polar_base[8]= 0.0;
+}
Modified: mc/1D/hc/trunk/hc_output.c
===================================================================
--- mc/1D/hc/trunk/hc_output.c 2011-07-16 19:19:23 UTC (rev 18778)
+++ mc/1D/hc/trunk/hc_output.c 2011-07-16 22:37:33 UTC (rev 18779)
@@ -387,3 +387,109 @@
}
fclose(out);
}
+/*
+
+print a simple VTK file given already expanded input
+(called from hc_print_spatial)
+
+
+*/
+
+void hc_print_vtk(FILE *out,float *xloc,float *xvec,
+ int npoints,int nlay,
+ hc_boolean binary)
+{
+ int i,ilay,ndata,poff;
+ hc_boolean little_endian;
+ /* determine machine type */
+ little_endian = hc_is_little_endian();
+ /* */
+ ndata = npoints*3;
+ /* print VTK */
+ fprintf(out,"# vtk DataFile Version 4.0\n");
+ fprintf(out,"generated by hc_print_vtk\n");
+ if(binary)
+ fprintf(out,"BINARY\n");
+ else
+ fprintf(out,"ASCII\n");
+ fprintf(out,"DATASET POLYDATA\n");
+
+ fprintf(out,"POINTS %i float\n",npoints*nlay);
+ for(ilay=0;ilay < nlay;ilay++){
+ for(i=0;i < npoints;i++){
+ poff = ilay * ndata + i*3;
+ if(binary)
+ 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]);
+ }
+ }
+ fprintf(out,"POINT_DATA %i\n",npoints*nlay);
+ fprintf(out,"VECTORS velocity float\n");
+ for(ilay=0;ilay < nlay;ilay++){
+ for(i=0;i < npoints;i++){
+ poff = ilay * ndata + i*3;
+ if(binary)
+ hc_print_be_float((xvec+poff),3,out,little_endian);
+ else
+ fprintf(out,"%.6e %.6e %.6e\n",xvec[poff],xvec[poff+1],
+ xvec[poff+2]);
+ }
+ }
+
+}
+/*
+ print big endian binary
+*/
+void hc_print_be_float(float *x, int n, FILE *out, hc_boolean little_endian)
+{
+ int i;
+ float *xcopy;
+ static size_t len = sizeof(float);
+ 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);
+ }else{
+ /* can write as is */
+ fwrite(x,len,n,out);
+ }
+}
+hc_boolean hc_is_little_endian(void)
+{
+ static const unsigned long a = 1;
+ return *(const unsigned char *)&a;
+}
+
+/*
+
+flip endianness of x
+
+*/
+void hc_flip_byte_order(void *x, size_t len)
+{
+ void *copy;
+ copy = (void *)malloc(len);
+ if(!copy){
+ fprintf(stderr,"flip_byte_order: memerror with len: %i\n",(int)len);
+ exit(-1);
+ }
+ memcpy(copy,x,len);
+ hc_flipit(x,copy,len);
+ free(copy);
+}
+/* this should not be called with (i,i,size i) */
+void hc_flipit(void *d, void *s, size_t len)
+{
+ unsigned char *dest = d;
+ unsigned char *src = s;
+ src += len - 1;
+ for (; len; len--)
+ *dest++ = *src--;
+}
+
Modified: mc/1D/hc/trunk/sh_exp.c
===================================================================
--- mc/1D/hc/trunk/sh_exp.c 2011-07-16 19:19:23 UTC (rev 18778)
+++ mc/1D/hc/trunk/sh_exp.c 2011-07-16 22:37:33 UTC (rev 18779)
@@ -1325,50 +1325,12 @@
float z, FILE *out)
{
int j,k;
- float xp[3],lon,lat;
+ float lon,lat;
for(j=0;j < exp[0].npoints;j++){
/*
get coordinates
*/
- switch(exp->type){
-#ifdef HC_USE_HEALPIX
-
- case SH_HEALPIX:
- switch(exp[0].heal.ordering){
- case SH_HEALPIX_RING:
- pix2ang_ring((long)exp[0].heal.nside,(long)j,(xp+HC_THETA),(xp+HC_PHI));
- break;
- case SH_HEALPIX_NEST:
- pix2ang_nest((long)exp[0].heal.nside,(long)j,(xp+HC_THETA),(xp+HC_PHI));
- break;
- default:
- fprintf(stderr,"print_sh_spatial_data: error: ordering %i undefined\n",
- exp[0].heal.ordering);
- exit(-1);
- break;
- }
- break; /* end Healpix branch */
-#endif
- case SH_RICK:
- /* compute location */
-#ifdef NO_RICK_FORTRAN
- rick_pix2ang(j,exp[0].lmax,(xp+HC_THETA),(xp+HC_PHI),
- &exp->rick);
-#else
- rick_f90_pix2ang(&j,&exp[0].lmax,(xp+HC_THETA),(xp+HC_PHI));
-#endif
- break;
-#ifdef HC_USE_SPHEREPACK
- case SH_SPHEREPACK_GAUSS:
- case SH_SPHEREPACK_EVEN:
- break;
-#endif
- default:
- sh_exp_type_error("sh_print_spatial_data",exp);
- break;
- } /* end type branch */
- lon = PHI2LON(xp[HC_PHI]);
- lat = THETA2LAT(xp[HC_THETA]);
+ sh_get_coordinates(exp,j,&lon,&lat);
/* print coordinates */
if(!use_3d){
/* print lon lat */
@@ -1383,6 +1345,53 @@
} /* end points in lateral space loop */
}
+void sh_get_coordinates(struct sh_lms *exp,
+ int i, float *lon, float *lat)
+{
+ float xp[3];
+ switch(exp->type){
+#ifdef HC_USE_HEALPIX
+
+ case SH_HEALPIX:
+ switch(exp[0].heal.ordering){
+ case SH_HEALPIX_RING:
+ pix2ang_ring((long)exp[0].heal.nside,
+ (long)i,(xp+HC_THETA),(xp+HC_PHI));
+ break;
+ case SH_HEALPIX_NEST:
+ pix2ang_nest((long)exp[0].heal.nside,
+ (long)i,(xp+HC_THETA),(xp+HC_PHI));
+ break;
+ default:
+ fprintf(stderr,"print_sh_spatial_data: error: ordering %i undefined\n",
+ exp[0].heal.ordering);
+ exit(-1);
+ break;
+ }
+ break; /* end Healpix branch */
+#endif
+ case SH_RICK:
+ /* compute location */
+#ifdef NO_RICK_FORTRAN
+ rick_pix2ang(i,exp[0].lmax,(xp+HC_THETA),(xp+HC_PHI),
+ &exp->rick);
+#else
+ rick_f90_pix2ang(&i,&exp[0].lmax,(xp+HC_THETA),(xp+HC_PHI));
+#endif
+ break;
+#ifdef HC_USE_SPHEREPACK
+ case SH_SPHEREPACK_GAUSS:
+ case SH_SPHEREPACK_EVEN:
+ break;
+#endif
+ default:
+ sh_exp_type_error("sh_print_spatial_data",exp);
+ break;
+ } /* end type branch */
+ *lon = PHI2LON(xp[HC_PHI]);
+ *lat = THETA2LAT(xp[HC_THETA]);
+}
+
/*
regular grid version
Modified: mc/1D/hc/trunk/sh_syn.c
===================================================================
--- mc/1D/hc/trunk/sh_syn.c 2011-07-16 19:19:23 UTC (rev 18778)
+++ mc/1D/hc/trunk/sh_syn.c 2011-07-16 22:37:33 UTC (rev 18779)
@@ -172,8 +172,7 @@
hc_svecalloc(&data,exp[0].npoints * shps,"sh_syn");
sh_compute_spatial(exp,ivec,FALSE,&dummy,data,verbose);
/* output */
- sh_print_spatial_data_to_file(exp,shps,data,
- (nset>1)?(TRUE):(FALSE),
+ sh_print_spatial_data_to_file(exp,shps,data,(nset>1)?(TRUE):(FALSE),
zlabel,stdout);
}
free(exp);free(data);
More information about the CIG-COMMITS
mailing list