[cig-commits] r11605 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Wed Mar 26 20:36:48 PDT 2008
Author: becker
Date: 2008-03-26 20:36:48 -0700 (Wed, 26 Mar 2008)
New Revision: 11605
Added:
mc/1D/hc/trunk/rotvec2vel.c
Modified:
mc/1D/hc/trunk/Makefile
mc/1D/hc/trunk/README.TXT
mc/1D/hc/trunk/ggrd_grdtrack_util.c
mc/1D/hc/trunk/ggrd_grdtrack_util.h
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/sh_ana.c
mc/1D/hc/trunk/sh_exp.c
Log:
Modified GMT4 handling, yet again. Added input from GMT netcdf grid.
Added rotvec2vel utility.
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2008-03-27 02:31:24 UTC (rev 11604)
+++ mc/1D/hc/trunk/Makefile 2008-03-27 03:36:48 UTC (rev 11605)
@@ -132,7 +132,7 @@
all: dirs libs hc hc_extract_sh_layer \
- sh_syn sh_ana sh_power
+ sh_syn sh_ana sh_power rotvec2vel
libs: dirs hc_lib $(HEAL_LIBS) $(RICK_LIB)
@@ -161,13 +161,17 @@
sh_power: $(LIBS) $(INCS) $(ODIR)/sh_power.o
$(CC) $(LIB_FLAGS) $(ODIR)/sh_power.o \
- -o $(BDIR)/sh_power -lhc -lrick $(HEAL_LIBS_LINKLINE) -lm $(LDFLAGS)
+ -o $(BDIR)/sh_power -lhc -lrick $(HEAL_LIBS_LINKLINE) $(GGRD_LIBS_LINKLINE) -lm $(LDFLAGS)
sh_ana: $(LIBS) $(INCS) $(ODIR)/sh_ana.o
$(CC) $(LIB_FLAGS) $(ODIR)/sh_ana.o \
- -o $(BDIR)/sh_ana -lhc -lrick $(HEAL_LIBS_LINKLINE) -lm $(LDFLAGS)
+ -o $(BDIR)/sh_ana -lhc -lrick $(HEAL_LIBS_LINKLINE) $(GGRD_LIBS_LINKLINE) -lm $(LDFLAGS)
+rotvec2vel: rotvec2vel.c
+ $(CC) $(CFLAGS) rotvec2vel.c -o $(BDIR)/rotvec2vel -lm $(LDFLAGS)
+
+
hc: $(LIBS) $(INCS) $(ODIR)/main.o $(PREM_OBJS)
$(CC) $(LIB_FLAGS) $(ODIR)/main.o -o $(BDIR)/hc \
-lhc -lrick $(HEAL_LIBS_LINKLINE) $(PREM_OBJS) \
Modified: mc/1D/hc/trunk/README.TXT
===================================================================
--- mc/1D/hc/trunk/README.TXT 2008-03-27 02:31:24 UTC (rev 11604)
+++ mc/1D/hc/trunk/README.TXT 2008-03-27 03:36:48 UTC (rev 11605)
@@ -219,7 +219,7 @@
harmonics expansion of degree 31
- a) obtain scalara data at the Gaussian intergration points
+ a) obtain scalar data at the Gaussian intergration points
sh_ana -127 | grdtrack -Lx -G$datadir/etopo5/etopo5.1.grd > etopo5.dat
@@ -228,6 +228,10 @@
cat etopo5.dat | sh_ana 127 > etopo5.ab
+ c) Alternative: use the grd file directly
+
+ sh_ana 127 $datadir/etopo5/etopo5.1.grd
+
2) convert spherical harmonics to spatial expansion
cat etopo5.ab | sh_syn > etopo5.127.dat
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2008-03-27 02:31:24 UTC (rev 11604)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2008-03-27 03:36:48 UTC (rev 11605)
@@ -390,7 +390,7 @@
char *grdfile, /* name, or prefix, of grd file with scalars */
struct GRD_HEADER **grd, /* pass as empty */
struct GMT_EDGEINFO **edgeinfo, /* pass as empty */
- char *edgeinfo_string, /* -L type flags from GMT, can be empty */
+ char *edgeinfo_string, /* -fg/ -L type flags from GMT, can be empty */
ggrd_boolean *geographic_in, /* this is normally TRUE */
int *pad, /* [4] array with padding (output) */
ggrd_boolean three_d, char *dfile, /* depth file name */
@@ -427,9 +427,6 @@
/*
init first edgeinfo
*/
-#ifdef USE_GMT4
- GMT_io_init (); /* Init the table i/o structure */
-#endif
GMT_boundcond_init (*edgeinfo);
/* check if geographic */
@@ -442,6 +439,11 @@
}else{
*geographic_in = FALSE;
}
+#ifdef USE_GMT4
+ GMT_io_init ();/* Init the table i/o structure */
+ GMT_grdio_init();
+ GMT_program = "g";
+#endif
if(verbose >= 2)
if(*geographic_in)
fprintf(stderr,"ggrd_grdtrack_init: detected geographic region from geostring: %s\n",
@@ -525,9 +527,10 @@
return 4;
}
#else /* 4.1.2 */
- sprintf((*grd)->name,"%s",grdfile);
- if (GMT_cdf_read_grd_info ((*grd))) {
- fprintf (stderr, "%s: error opening file %s\n",
+ if(GMT_read_grd_info (grdfile,*grd)){
+
+ //if (GMT_cdf_read_grd_info ((*grd))) {
+ fprintf (stderr, "%s: error opening file %s for header\n",
"ggrd_grdtrack_init", grdfile);
return 4;
}
@@ -542,9 +545,8 @@
"ggrd_grdtrack_init", filename);
return 6;
}
-#else
- sprintf((*grd+i)->name,"%s",filename);
- if (GMT_cdf_read_grd_info ((*grd+i))) {
+#else /* gmt 4 */
+ if (GMT_read_grd_info (filename,(*grd+i))) {
fprintf (stderr, "%s: error opening file %s (-D option was used)\n",
"ggrd_grdtrack_init", filename);
return 6;
@@ -643,14 +645,15 @@
read the grd files
*/
#ifdef USE_GMT4
- sprintf((*grd+i)->name,"%s",filename);
- if (GMT_cdf_read_grd ((*grd+i), (*f+i* (*mm)),
+ /* GMT 4 */
+ if (GMT_read_grd (filename,(*grd+i), (*f+i* (*mm)),
*west, *east, *south, *north,
pad, FALSE)) {
fprintf (stderr, "%s: error reading file %s\n", "ggrd_grdtrack_init", grdfile);
return 10;
}
#else
+ /* old GMT */
if (GMT_cdf_read_grd (filename, (*grd+i), (*f+i* (*mm)),
*west, *east, *south, *north,
pad, FALSE)) {
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.h
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.h 2008-03-27 02:31:24 UTC (rev 11604)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.h 2008-03-27 03:36:48 UTC (rev 11605)
@@ -8,6 +8,7 @@
#ifndef __GGRD_READ_GMT__
#include "gmt.h"
+void GMT_grdio_init (void);
#define __GGRD_READ_GMT__
#endif
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2008-03-27 02:31:24 UTC (rev 11604)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2008-03-27 03:36:48 UTC (rev 11605)
@@ -1,4 +1,5 @@
/* ggrd_grdtrack_util.c */
+void ggrd_init_master(struct ggrd_master *);
int ggrd_grdtrack_init_general(unsigned char, char *, char *, char *, struct ggrd_gt *, unsigned char, unsigned char);
int ggrd_grdtrack_rescale(struct ggrd_gt *, unsigned char, unsigned char, unsigned char, double);
unsigned char ggrd_grdtrack_interpolate_rtp(double, double, double, struct ggrd_gt *, double *, unsigned char);
@@ -12,8 +13,11 @@
void ggrd_gt_interpolate_z(double, float *, int, int *, int *, double *, double *, unsigned char, unsigned char *);
void ggrd_interpol_time(double, struct ggrd_t *, int *, int *, double *, double *, double);
int interpolate_seafloor_ages(double, double, double, struct ggrd_master *, double *);
+FILE *ggrd_open(char *, char *, char *);
void ggrd_vecalloc(double **, int, char *);
void ggrd_vecrealloc(double **, int, char *);
+void ggrd_calc_mean_and_stddev(double *, double *, int, double *, double *, double *, unsigned char, unsigned char, double *);
+void ggrd_indexx(int, double *, int *);
float ggrd_gt_rms(float *, int);
float ggrd_gt_mean(float *, int);
/* ggrd_readgrds.c */
@@ -49,7 +53,6 @@
void hc_ludcmp_3x3(double [3][3], int *);
void hc_lubksb_3x3(double [3][3], int *, double *);
/* hc_misc.c */
-FILE *hc_open(char *, char *, char *);
void hc_dvecalloc(double **, int, char *);
void hc_svecalloc(float **, int, char *);
void hc_vecalloc(double **, int, char *);
@@ -69,8 +72,6 @@
char *hc_name_boolean(unsigned short);
unsigned short hc_toggle_boolean(unsigned short *);
void hc_advance_argument(int *, int, char **);
-void hc_calc_mean_and_stddev(double *, double *, int, double *, double *, double *, unsigned short, unsigned short, double *);
-void hc_indexx(int, double *, int *);
/* 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);
@@ -151,6 +152,8 @@
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);
Added: mc/1D/hc/trunk/rotvec2vel.c
===================================================================
--- mc/1D/hc/trunk/rotvec2vel.c (rev 0)
+++ mc/1D/hc/trunk/rotvec2vel.c 2008-03-27 03:36:48 UTC (rev 11605)
@@ -0,0 +1,286 @@
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+/*
+
+utility program, here for historical reasons
+
+
+ */
+#define DEF_FIXED -1
+
+//
+// reads in rotation vector file and xy code file with plate codes
+//
+// output is lon lat vp vt
+//
+// these are the coded plates as in the lon lat code tripel that is
+// read from stdin
+//
+#define CODED_PLATES 14
+#define PLATE_CODES "ANT AUS AFR PAC EUR NAM NAZ COC CAR ARA PHI SAM IND JDF "
+
+// some factors
+#define PIOVERONEEIGHTY 0.0174532925199433
+#define PIHALF 1.5707963267949
+#define PI 3.14159265358979
+//#define VELFACTOR (6371009.0/1.0e6*1.0e2)
+//#define PIOVERONEEIGHTY_TIMES_VELFACTOR (PIOVERONEEIGHTY*VELFACTOR)
+#define PIOVERONEEIGHTY_TIMES_VELFACTOR 11.1195083724191
+FILE *myopen(const char *,const char *);
+
+
+int main(int argc, char *argv[])
+{
+ int i,j,k,nplt,nrp,fixed_plate,*stats,minplate,maxplate,normalize=0,unity_vec=0,
+ assigned[CODED_PLATES],hit,name[CODED_PLATES],code,code3;
+ double *rotvel,*points,x,y,z,stheta,cphi,ctheta,sphi,vx,vy,vz,vt,vp,
+ redvec[3],length;
+ char allplates[]=PLATE_CODES,pname[4],tmpstring[300];
+ FILE *vel;
+
+ switch(argc){
+ case 2:{
+ fixed_plate=DEF_FIXED;
+ break;
+ }
+ case 3:{
+ sscanf(argv[2],"%i",&fixed_plate);
+ break;
+ }
+ default:{
+ fprintf(stderr,"%s rotvector_file [fixed_plate]\n",argv[0]);
+ fprintf(stderr,"\t reads points in (x y plate_code) format from stdin\n");
+ fprintf(stderr,"\t and calculates velocities (vp,vt) from rotvector_file\n");
+ fprintf(stderr,"\t rotvectorfile is in\n");
+ fprintf(stderr,"\t platename_1 wx_1 wy_1 wz_1\n");
+ fprintf(stderr,"\t platename_2 wx_2 wy_2 wz_2...\n");
+ fprintf(stderr,"\t format, wi_j in input is in [degrees/Myr]\n");
+ fprintf(stderr,"\t fixed_plate is the reference plate number (1,2,...) in the poles list,\n");
+ fprintf(stderr,"\t if set to -1, no change in motions.\n");
+ fprintf(stderr,"\t if set to -2, normalize all omega vectors to 1 deg/Myr length\n");
+ fprintf(stderr,"\t if set to -3, use 1 deg/Myr w_x vector, rest zero\n");
+ fprintf(stderr,"\t if set to -4, use 1 deg/Myr w_y vector, rest zero\n");
+ fprintf(stderr,"\t if set to -5, use 1 deg/Myr w_z vector, rest zero\n");
+ fprintf(stderr,"\t By default set to %i.\n",DEF_FIXED);
+ fprintf(stderr,"\t output is in\n\tlon lat v_p v_t\n\tformat in cm/yr\n");
+ exit(-1);
+ break;
+ }}
+ for(i=0;i<CODED_PLATES;i++)
+ assigned[i]=0;
+
+ vel=myopen(argv[1],"r");
+
+ rotvel=(double *)malloc(sizeof(double)*3*CODED_PLATES);
+ stats=(int *)calloc(CODED_PLATES,sizeof(int));
+
+ fprintf(stderr,"%s: reading from %s, fixed plate: %i\n",argv[0],argv[1],fixed_plate);
+ nplt=0;
+ // skip initial comment lines
+ while(fscanf(vel,"%[^\n]%*1c",tmpstring)!=EOF){
+ if(tmpstring[0]!='#' && tmpstring[0]!=' '){
+ if(sscanf(tmpstring,"%s %lf %lf %lf",pname,&vx,&vy,&vz)==4){
+ for(hit=i=0;i<CODED_PLATES;i++)
+ if(strncmp(pname,(allplates+4*i),3)==0){
+ if(assigned[i]){
+ fprintf(stderr,"%s: plate %s was assigned already, is it twice in list?\n",
+ argv[0],pname);
+ exit(-1);
+ }
+ assigned[i]=hit=1;
+ k=i*3;
+ *(rotvel+k++)=vx;
+ *(rotvel+k++)=vy;
+ *(rotvel+k)=vz;
+
+ name[nplt]=i;
+ }
+ if(!hit){
+ fprintf(stderr,"%s: could not find internal code for input plate code %s\n",
+ argv[0],pname);
+ fprintf(stderr,"%s: will ignore this code\n",argv[0]);
+ fprintf(stderr,"%s: internal codes are:\n",argv[0]);
+ for(j=0;j<CODED_PLATES;j++){
+ strncpy(pname,(allplates+j*4),3);
+ fprintf(stderr,"%s: plate nr %i, code: %s\n",argv[0],j+1,pname);
+ }
+ }
+ nplt++;
+ }
+ }
+ }
+ fclose(vel);
+ fprintf(stderr,"%s: read in %i poles\n",argv[0],nplt);
+
+ // assign zero rotation vectors to all plates we did not find in the list
+ for(i=0,k=0;i<CODED_PLATES;i++,k+=3)
+ if(!assigned[i])
+ for(j=0;j<3;j++)
+ rotvel[k+j]=0.0;
+
+ if(fixed_plate>0){
+ fprintf(stderr,"%s: input plate number %i will be fixed\n",argv[0],fixed_plate);
+ if(fixed_plate>nplt){
+ fprintf(stderr,"%s: fixed plate number (%i) greater than plate numbers (1...%i)\n\n",
+ argv[0],fixed_plate,nplt);
+ }
+ for(i=0;i<3;i++)
+ redvec[i]= *(rotvel+name[fixed_plate-1]*3+i);
+ }else{
+ for(i=0;i<3;i++)
+ redvec[i]=0.0;
+ if(fixed_plate==-2){
+ normalize=1;
+ fprintf(stderr,"%s: normalizing all vectors\n",argv[0]);
+ }else if(fixed_plate < -2 && fixed_plate > -6){
+ unity_vec=-fixed_plate-2;
+ fprintf(stderr,"%s: using unity for %ith component of omega, zero else\n",
+ argv[0],unity_vec);
+ }
+ }
+ //
+ // convert rotation poles
+ //
+ for(k=i=0;i<CODED_PLATES;i++,k+=3){
+ if(assigned[i]){
+ for(j=0;j<3;j++)
+ *(rotvel+k+j) -= redvec[j];
+ strncpy(pname,(allplates+4*i),3);
+ // original vector
+ fprintf(stderr,"%s: plate %3i %3s wx: %12g wy: %12g wz: %12g",
+ argv[0],i+1,pname,*(rotvel+k),*(rotvel+k+1),*(rotvel+k+2));
+ if(i == name[fixed_plate-1])
+ fprintf(stderr," [deg/Myr] (fixed)\n");
+ else
+ fprintf(stderr," [deg/Myr]\n");
+
+ // normalize
+ for(length=0.0,j=0;j<3;j++)
+ length+= *(rotvel+k+j)* *(rotvel+k+j);
+ length=sqrt(length);
+ if(normalize)
+ for(j=0;j<3;j++)
+ *(rotvel+k+j) /= length;
+ // or assign unit length
+ if(unity_vec){
+ for(j=0;j<3;j++){
+ if(unity_vec-1==j){
+ *(rotvel+k+j)=1.0;
+ }else{
+ *(rotvel+k+j)=0.0;
+ }
+ }
+ }
+ for(length=0.0,j=0;j<3;j++)
+ length+= *(rotvel+k+j)* *(rotvel+k+j);
+ length=sqrt(length);
+ // converted
+ fprintf(stderr,"%s: wx: %12g wy: %12g wz: %12g [deg/Myr]\n",
+ argv[0],*(rotvel+k),*(rotvel+k+1),*(rotvel+k+2));
+ fprintf(stderr,"%s: length: %12g [deg/Myr]\n",argv[0],length);
+
+ // rescale the omega vector such that it will yield cm/yr when
+ // multiplied with r vector of unit length
+ *(rotvel+k) *=PIOVERONEEIGHTY_TIMES_VELFACTOR;
+ *(rotvel+k+1)*=PIOVERONEEIGHTY_TIMES_VELFACTOR;
+ *(rotvel+k+2)*=PIOVERONEEIGHTY_TIMES_VELFACTOR;
+ }
+ }
+ minplate=1e6;
+ maxplate=-1;
+ // read in lon lat code tripels
+ nrp=1;
+ points=(double *)malloc(sizeof(double)*3*nrp);
+ k=0;
+ while(fscanf(stdin,"%lf %lf %lf",(points+k),(points+k+1),(points+k+2))==3){
+ if(*(points+k)<0.0)
+ *(points+k) += 360.0;
+ *(points+k) *=PIOVERONEEIGHTY;
+ *(points+k+1)*=PIOVERONEEIGHTY;
+ /* fix poles */
+ *(points+k+1)= PIHALF- *(points+k+1);
+ if(fabs(points[k+1]) < 1e-5)
+ points[k+1] = 1e-5;
+ if(fabs(points[k+1] -PIHALF) < 1e-5)
+ points[k+1] = PIHALF - 1e-5;
+
+ // change code here fore more efficient assignment
+ *(points+k+2)-=1;
+ if(*(points+k+2)<minplate)
+ minplate= *(points+k+2);
+ if(*(points+k+2)>maxplate)
+ maxplate= *(points+k+2);
+ nrp++;
+ k += 3;
+ if((points=((double *)realloc(points,sizeof(double)*3*nrp)))==NULL){
+ fprintf(stderr,"%s: memerror, too many points for x array, %i\n",argv[0],nrp);
+ exit(-1);
+ }
+ }
+ nrp--;
+
+ fprintf(stderr,"%s: read %i points for velocities, calculating...\n",argv[0],nrp);
+ fprintf(stderr,"%s: input plate codes run from %3i to %3i\n",argv[0],minplate+1,maxplate+1);
+ if(minplate<0){
+ fprintf(stderr,"%s: expect code to start from 1\n",argv[0]);exit(-1);
+ }
+ if(maxplate>nplt-1 || maxplate-minplate>nplt){
+ fprintf(stderr,"%s: could only read %i plate codes from rotation vector file\n",
+ argv[0],nplt);
+ fprintf(stderr,"%s: will set velocities for all other codes to zero\n",argv[0]);
+ }
+
+ fprintf(stderr,"%s: now output...\n",argv[0]);
+
+#define WX ( *(rotvel+code3) )
+#define WY ( *(rotvel+code3+1) )
+#define WZ ( *(rotvel+code3+2) )
+
+ for(i=k=0;i<nrp;i++,k+=3){
+ code=(int)*(points+k+2);
+ if(assigned[code]){
+ code3 = code * 3;
+ // we have included the radius factor in the omega vector,
+ // hence all locations are calculated with radius unity
+ x=(cphi=cos(*(points+k)))*(stheta=sin(*(points+k+1)));
+ y=(sphi=sin(*(points+k)))*stheta;
+ z=(ctheta=cos(*(points+k+1)));
+ vx = WY * z - WZ * y;
+ vy = WZ * x - WX * z;
+ vz = WX * y - WY * x;
+ *(stats+code)+=1;
+ /* v_theta goes southward */
+ vt= ctheta*cphi*vx + ctheta*sphi*vy - stheta*vz;
+ /* v_phi goes eastward */
+ vp= -sphi*vx + cphi*vy;
+ fprintf(stdout,"%12g %12g %20g %20g\n",
+ (1.0/PIOVERONEEIGHTY)* *(points+k),
+ 90.0-(1.0/PIOVERONEEIGHTY)* *(points+k+1),
+ vp,vt);
+ }else{
+ ;
+ //fprintf(stdout,"%12g %12g %20g %20g\n",
+ // (1.0/PIOVERONEEIGHTY)* *(points+k),
+ // 90.0-(1.0/PIOVERONEEIGHTY)* *(points+k+1),0.,0.);
+ }
+ }
+ for(i=0;i<CODED_PLATES;i++)
+ if(assigned[i]){
+ strncpy(pname,(allplates+4*i),3);
+ fprintf(stderr,"%s: %3s (%3i) was assigned %7i times out of %7i, %15.8f percent\n",
+ argv[0],pname,i+1,*(stats+i),nrp,(double)(*(stats+i))/(double)nrp*100.0);
+ }
+ fprintf(stderr,"%s: done\n",argv[0]);
+ return 0;
+}
+FILE *myopen(const char *name,const char *modus)
+{
+ FILE *tmp;
+ if( (tmp = (FILE *)fopen(name,modus)) == NULL)
+ { fprintf(stderr,"Error opening \"%s\" for \"%s\".\n",name,modus);
+ fprintf(stderr,"Exiting.\n");
+ exit(-1);};
+ return ((FILE *)tmp);
+}
Modified: mc/1D/hc/trunk/sh_ana.c
===================================================================
--- mc/1D/hc/trunk/sh_ana.c 2008-03-27 02:31:24 UTC (rev 11604)
+++ mc/1D/hc/trunk/sh_ana.c 2008-03-27 03:36:48 UTC (rev 11605)
@@ -10,14 +10,27 @@
Thorsten Becker (twb at usc.edu)
-usage: cat data.lonlatz | sh_ana l_max ivec
+usage:
+ cat data.lonlatz | sh_ana l_max ivec type short_format
+
l_max: max order of expansion. if negative, will print out the spatial
locations needed on input
ivec: 0: expand scalar field
1: expand vector field
+ file: if filename other than 0, 1, or vec_t.grd will read from Netcdf/GMT grid file
+ vec_t.grd: will read theta and phi components of vector field from vec_t.grd and vec_p.grd
+ type: 0: Rick's spherical harmonics
+ 1: Healpix (defunct)
+
+ short_format: header format
+
+ 0: long output for hc
+ 1: short output, only lmax
+
+
where data.lonlatz has the data at the required locations as put out by sh_ana -lmax
@@ -38,13 +51,15 @@
int main(int argc, char **argv)
{
- int type = SH_RICK,lmax,shps, nset=1,ivec=0,ilayer=0,i;
+ int type = SH_RICK,lmax,shps, nset=1,ivec=0,ilayer=0,i,failed;
struct sh_lms *exp;
- hc_boolean verbose = TRUE, use_3d = FALSE, short_format = FALSE,
+ 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;
+ struct ggrd_gt ggrd[2];
SH_RICK_PREC *dummy;
HC_PREC fac[3] = {1.,1.,1.};
+ char cdummy;
/*
command line parameters
@@ -60,8 +75,37 @@
}
}
}
- if(argc > 2)
- sscanf(argv[2],"%i",&ivec);
+ if(argc > 2){
+ if((strcmp(argv[2],"0")==0)||(strcmp(argv[2],"1")==0)) /* regular input */
+ sscanf(argv[2],"%i",&ivec);
+ else{
+ /* GMT grd file input */
+#ifdef USE_GMT4
+ failed = ggrd_grdtrack_init_general(FALSE,argv[2],&cdummy,"-fg",ggrd,TRUE,FALSE);
+#else
+ failed = ggrd_grdtrack_init_general(FALSE,argv[2],&cdummy,"-Lg",ggrd,TRUE,FALSE);
+#endif
+ if(failed){
+ fprintf(stderr,"%s: error opening netcdf grd file\n",argv[0]);
+ exit(-1);
+ }
+ if(strcmp(argv[2],"vec_t.grd")==0){ /* vectors */
+#ifdef USE_GMT4
+ failed = ggrd_grdtrack_init_general(FALSE,"vec_p.grd",&cdummy,"-fg",(ggrd+1),TRUE,FALSE);
+#else
+ failed = ggrd_grdtrack_init_general(FALSE,"vec_p.grd",&cdummy,"-Lg",(ggrd+1),TRUE,FALSE);
+#endif
+ if(failed){
+ fprintf(stderr,"%s: error opening second netcdf grd file %s\n",argv[0],"vec_p.grd");
+ exit(-1);
+ }
+ ivec = 1;
+ }else{
+ ivec = 0;
+ }
+ read_grd = TRUE;
+ }
+ }
if(argc > 3)
sscanf(argv[3],"%i",&type);
if(argc > 4){
@@ -69,13 +113,15 @@
short_format = (hc_boolean)i;
}
if((argc > 4)||(argc<=1)){
- fprintf(stderr,"usage: %s l_max [ivec, %i] [type, %i] [short_format, %i\n",
+ fprintf(stderr,"usage: %s l_max [ivec, %i] [type, %i] [short_format, %i]\n",
argv[0],ivec,type,short_format);
fprintf(stderr," l_max: max order of expansion. if negative, will print out the spatial\n");
fprintf(stderr," locations needed on input\n");
fprintf(stderr," for Rick SH format, lmax needs to be 2**n-1\n\n");
fprintf(stderr," ivec: 0: expand scalar field (input: lon lat scalar)\n");
- fprintf(stderr," 1: expand vector field (input: lon lat v_t v_p\n\n");
+ fprintf(stderr," 1: expand vector field (input: lon lat v_t v_p\n");
+ fprintf(stderr," vec_t.grd: expand vector field in vec_t.grd (theta) and vec_p.grd (phi)\n");
+ fprintf(stderr," file.grd: expand scalar in file.grd, where .grd indicates global Netcdf files\n\n");
fprintf(stderr," type %i: use Rick's routines internally (output format DT convection)\n",SH_RICK);
fprintf(stderr," %i: use Healpix's routines internally (output format DT convection)\n\n",SH_HEALPIX);
fprintf(stderr,"short_format: 0: use long header files format as in HC\n");
@@ -85,9 +131,16 @@
if(print_spatial_base)
fprintf(stderr,"%s: printing spatial base for lmax: %i\n",
argv[0],lmax);
- else
- fprintf(stderr,"%s: expanding to lmax: %i, expecting %s\n",
- argv[0],lmax,(ivec)?("lon lat vt vp"):("lon lat x"));
+ else{
+ if(read_grd)
+ fprintf(stderr,"%s: expanding to lmax: %i, reading from %s %s\n",
+ argv[0],lmax,argv[2],(ivec)?("vec_p.grd"):(""));
+
+ else
+ fprintf(stderr,"%s: expanding to lmax: %i, expecting %s\n",
+ argv[0],lmax,(ivec)?("lon lat vt vp"):("lon lat x"));
+
+ }
/*
select numbers of expansions, scalar or pol/tor for vector field
*/
@@ -105,8 +158,12 @@
/*
perform spherical harmonic analysis
*/
- /* read in data from stdin */
- sh_read_spatial_data_from_file(exp,stdin,use_3d,shps,data,&zlabel);
+ if(read_grd){
+ sh_read_spatial_data_from_grd(exp,ggrd,use_3d,shps,data,&zlabel);
+ }else{
+ /* read in data from stdin */
+ sh_read_spatial_data_from_file(exp,stdin,use_3d,shps,data,&zlabel);
+ }
/*
perform spherical harmonic expansion
*/
Modified: mc/1D/hc/trunk/sh_exp.c
===================================================================
--- mc/1D/hc/trunk/sh_exp.c 2008-03-27 02:31:24 UTC (rev 11604)
+++ mc/1D/hc/trunk/sh_exp.c 2008-03-27 03:36:48 UTC (rev 11605)
@@ -664,8 +664,29 @@
my_boolean use_3d,
int shps, float *data, float *z)
{
+ struct ggrd_gt *gdummy=(struct ggrd_gt *)NULL;
+ sh_read_spatial_data(exp,in,gdummy,FALSE,use_3d,shps,data,z);
+}
+/* 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)
+{
+ FILE *fdummy=NULL;
+ sh_read_spatial_data(exp,fdummy,ggrd,TRUE,use_3d,shps,data,z);
+}
+/*
+
+generic function
+
+*/
+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)
+
+{
float lon,lat,xp[3];
int j,k;
+ double dvalue;
/*
read in data for each layer
*/
@@ -715,51 +736,71 @@
read coordinates
*/
if(!use_3d){
+ if(!use_grd){
+ /*
+ read in lon lat
+ */
+ if(fscanf(in,"%f %f",&lon,&lat) != 2){
+ fprintf(stderr,"sh_read_spatial_data: error: lon lat format: pixel %i: read error\n",
+ (int)j);
+ exit(-1);
+ }
+ }
+ }else{
+ if(!use_grd){
+ /*
+ read in lon lat z[i]
+ */
+ if(fscanf(in,"%f %f %f",&lon,&lat,z) != 3){
+ fprintf(stderr,"sh_read_spatial_data: error: lon lat z format: pixel %i: read error\n",
+ (int)j);
+ exit(-1);
+ }
+ }else{
+ fprintf(stderr,"sh_read_spatial_data: error: 3D not implemented for grd yet\n");
+ exit(-1);
+ }
+ }
+ if(use_grd){
/*
- read in lon lat
+ interpolate from grd
*/
- if(fscanf(in,"%f %f",&lon,&lat) != 2){
- fprintf(stderr,"sh_read_spatial_data: error: lon lat format: pixel %i: read error\n",
- (int)j);
- exit(-1);
+ for(k=0;k < shps;k++){
+ if(!ggrd_grdtrack_interpolate_tp((double)xp[HC_THETA],(double)xp[HC_PHI],(ggrd+k),&dvalue,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]));
+ exit(-1);
+ }
+ data[k*exp[0].npoints+j] = dvalue;
}
}else{
/*
- read in lon lat z[i]
+ read data
*/
- if(fscanf(in,"%f %f %f",&lon,&lat,z) != 3){
- fprintf(stderr,"sh_read_spatial_data: error: lon lat z format: pixel %i: read error\n",
- (int)j);
- exit(-1);
+ for(k=0;k < shps;k++){
+ if(fscanf(in,"%f",(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);
+ }
}
- }
- /*
- read data
- */
- for(k=0;k < shps;k++){
- if(fscanf(in,"%f",(data+k*exp[0].npoints+j))!=1){
- fprintf(stderr,"sh_read_spatial_data: error: scalar format: pixel %i: read error\n",
+ /*
+ adjust longitude range
+ */
+ if(lon < 0)
+ lon += 360.0;
+ /*
+ check if location is OK (we don't know about z defaults)
+ */
+ if(((fabs(PHI2LON(xp[HC_PHI])-lon) > 1e-3)&&(fabs(fabs(PHI2LON(xp[HC_PHI])-lon)-360) > 1e-3))||
+ (fabs(THETA2LAT(xp[HC_THETA])-lat) > 1e-3)){
+ 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);
exit(-1);
}
}
- /*
- adjust longitude range
- */
- if(lon < 0)
- lon += 360.0;
- /*
- check if location is OK (we don't know about z defaults)
- */
- if(((fabs(PHI2LON(xp[HC_PHI])-lon) > 1e-3)&&
- (fabs(fabs(PHI2LON(xp[HC_PHI])-lon)-360) > 1e-3))||
- (fabs(THETA2LAT(xp[HC_THETA])-lat) > 1e-3)){
- 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);
- exit(-1);
- }
} /* end points in layer loop */
}
More information about the cig-commits
mailing list