[cig-commits] r11651 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Sun Mar 30 21:36:28 PDT 2008
Author: becker
Date: 2008-03-30 21:36:28 -0700 (Sun, 30 Mar 2008)
New Revision: 11651
Modified:
mc/1D/hc/trunk/ggrd_grdtrack_util.c
mc/1D/hc/trunk/ggrd_grdtrack_util.h
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/make_tar
Log:
Minor improvements for GMT4 handling.
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2008-03-31 04:21:42 UTC (rev 11650)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2008-03-31 04:36:28 UTC (rev 11651)
@@ -44,6 +44,18 @@
+/* for debugging, mostly */
+
+void ggrd_grdinfo(char *filename)
+{
+ struct ggrd_gt g[1];
+ char cdummy='c';
+ ggrd_grdtrack_init_general(FALSE,filename,&cdummy,"",g,2,FALSE);
+ fprintf(stderr,"ggrd_grdinfo: %s W: %g E: %g S: %g N: %g\n",
+ filename,g->grd->x_min,g->grd->x_max,g->grd->y_min,g->grd->y_max);
+
+}
+
/*
init structure and files for the general case, this is a wrapper
for what's below
@@ -418,34 +430,42 @@
float dz1,dz2;
struct GRD_HEADER ogrd;
int i,one_or_zero,nx,ny,mx,my,nn;
- char filename[BUFSIZ*2];
+ char filename[BUFSIZ*2],**cdummy;
+ static int gmt_init = FALSE;
/*
deal with edgeinfo
*/
*edgeinfo = (struct GMT_EDGEINFO *)
GMT_memory (VNULL, (size_t)1, sizeof(struct GMT_EDGEINFO), "ggrd_grdtrack_init");
+
+#ifdef USE_GMT4
+
+ if(!gmt_init){
+ /* this should be OK as is. init only once globally */
+ GMT_program = "ggrd";
+ GMT_make_fnan (GMT_f_NaN);
+ GMT_make_dnan (GMT_d_NaN);
+ GMT_io_init ();/* Init the table i/o structure */
+ GMT_grdio_init();
+ gmt_init = TRUE;
+ }
+
+#endif
/*
- init first edgeinfo
+ init first edgeinfo (period/global?)
*/
GMT_boundcond_init (*edgeinfo);
-
/* check if geographic */
- if (strlen(edgeinfo_string)>2){
+ if (strlen(edgeinfo_string)>2){ /* the boundary flag was set */
+ /* parse */
GMT_boundcond_parse (*edgeinfo, (edgeinfo_string+2));
- if ((*edgeinfo+0)->gn)
+ if ((*edgeinfo)->gn)
*geographic_in = TRUE;
else
*geographic_in = FALSE;
}else{
*geographic_in = FALSE;
}
-#ifdef USE_GMT4
- GMT_io_init ();/* Init the table i/o structure */
- GMT_grdio_init();
- GMT_program = "g";
- GMT_make_fnan (GMT_f_NaN);
- GMT_make_dnan (GMT_d_NaN);
-#endif
if(verbose >= 2)
if(*geographic_in)
fprintf(stderr,"ggrd_grdtrack_init: detected geographic region from geostring: %s\n",
@@ -519,22 +539,30 @@
GMT_memory (*edgeinfo, (size_t)(*nz), sizeof(struct GMT_EDGEINFO), "ggrd_grdtrack_init");
if(verbose >= 2)
fprintf(stderr,"ggrd_grdtrack_init: mem alloc ok\n");
+#ifdef USE_GMT4
+ /* init the header */
+ GMT_grd_init (*grd,0,cdummy,FALSE);
+#endif
if(*nz == 1){
if(verbose >= 2)
- fprintf(stderr,"ggrd_grdtrack_init: opening single file %s\n",grdfile);
+
#ifndef USE_GMT4 /* old */
+ fprintf(stderr,"ggrd_grdtrack_init: opening single file %s, GMT3 mode\n",grdfile);
if (GMT_cdf_read_grd_info (grdfile,(*grd))) {
fprintf (stderr, "%s: error opening file %s\n",
"ggrd_grdtrack_init", grdfile);
return 4;
}
+
#else /* 4.1.2 */
+ if(verbose >= 2)
+ fprintf(stderr,"ggrd_grdtrack_init: opening single file %s, GMT4 mode\n",grdfile);
if(GMT_read_grd_info (grdfile,*grd)){
fprintf (stderr, "%s: error opening file %s for header\n",
"ggrd_grdtrack_init", grdfile);
return 4;
}
-#endif
+#endif
}else{
/* loop through headers for testing purposess */
for(i=0;i<(*nz);i++){
@@ -625,10 +653,6 @@
pad on sides
*/
pad[0] = pad[1] = pad[2] = pad[3] = 2;
- if (verbose)
- fprintf(stderr,"ggrd_grdtrack_init: reading grd files (%g - %g (%i) %g - %g (%i); %i %s)\n",
- *west,*east,nx,*south,*north,ny,
- *geographic_in,(edgeinfo_string+2));
for(i=0;i < (*nz);i++){
/*
loop through layers
@@ -640,6 +664,10 @@
}else{ /* construct full filename */
sprintf(filename,"%s.%i.grd",grdfile,i+1);
}
+ if (verbose)
+ fprintf(stderr,"ggrd_grdtrack_init: reading grd file %s (%g - %g (%i) %g - %g (%i); geo: %i flag: %s\n",
+ filename,*west,*east,nx,*south,*north,ny,
+ *geographic_in,edgeinfo_string);
/*
read the grd files
@@ -647,11 +675,12 @@
#ifdef USE_GMT4
/* GMT 4 */
if (GMT_read_grd (filename,(*grd+i), (*f+i* (*mm)),
- *west, *east, *south, *north,
- pad, FALSE)) {
+ *west, *east, *south, *north,
+ pad, FALSE)) {
fprintf (stderr, "%s: error reading file %s\n", "ggrd_grdtrack_init", grdfile);
return 10;
}
+ //fprintf(stderr,"%g %g %i %i %i %i\n",(*grd)->z_scale_factor,(*grd)->z_add_offset,nx,ny,mx,my);
#else
/* old GMT */
if (GMT_cdf_read_grd (filename, (*grd+i), (*f+i* (*mm)),
@@ -679,21 +708,39 @@
#endif
}
/* Set boundary conditions */
- GMT_boundcond_set ((*grd+i), (*edgeinfo+i), pad,
- (*f+i*(*mm)));
+ GMT_boundcond_set ((*grd+i), (*edgeinfo+i), pad, (*f+i*(*mm)));
} /* end layer loop */
if(verbose){
- ggrd_print_layer_avg(*f,*z,*mm,*nz,stderr);
+ ggrd_print_layer_avg(*f,*z,mx,my,*nz,stderr,pad);
}
return 0;
}
-void ggrd_print_layer_avg(float *x,float *z,int n, int m,FILE *out)
+void ggrd_print_layer_avg(float *x,float *z,int nx, int ny, int m,FILE *out,int *pad)
{
- int i;
- for(i=0;i < m;i++){
- fprintf(stderr,"ggrd_grdtrack_init: layer %3i at depth %11g, mean: %11g rms: %11g\n",
- i+1,z[i],ggrd_gt_mean((x+i*n),n),ggrd_gt_rms((x+i*n),n));
+ int i,j,k,yl,xl,l,nxny,nxnyr;
+ float *tmp;
+ nxny = nx*ny; /* size with padding */
+ if(pad[0]+pad[1]+pad[2]+pad[3] == 0){
+ for(i=0;i < m;i++){
+ fprintf(stderr,"ggrd_grdtrack_init: layer %3i at depth %11g, mean: %11g rms: %11g\n",
+ i+1,z[i],ggrd_gt_mean((x+i*nxny),nxny),ggrd_gt_rms((x+i*nxny),nxny));
+ }
+ }else{
+
+ nxnyr = (nx-pad[0]-pad[1]) * (ny-pad[2]-pad[3]); /* actual data */
+ tmp = (float *)malloc(sizeof(float)*nxnyr);
+ if(!tmp)GGRD_MEMERROR("ggrd_print_layer_avg");
+ xl = nx - pad[1];
+ yl = ny - pad[2];
+ for(i=0;i < m;i++){ /* loop through depths */
+ for(l=0,j=pad[3];j < yl;j++)
+ for(k=pad[0];k < xl;k++,l++)
+ tmp[l] = x[i*nxny + j * nx + k];
+ fprintf(stderr,"ggrd_grdtrack_init: layer %3i at depth %11g, mean: %11g rms: %11g\n",
+ i+1,z[i],ggrd_gt_mean(tmp,nxnyr),ggrd_gt_rms(tmp,nxnyr));
+ }
+ free(tmp);
}
}
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.h
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.h 2008-03-31 04:21:42 UTC (rev 11650)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.h 2008-03-31 04:36:28 UTC (rev 11651)
@@ -120,7 +120,7 @@
float ggrd_gt_mean(float *,int );
FILE *ggrd_open(char *, char *, char *);
-void ggrd_print_layer_avg(float *,float *,int , int ,FILE *);
+void ggrd_print_layer_avg(float *,float *,int , int ,int, FILE *,int *);
void ggrd_find_spherical_vel_from_rigid_cart_rot(double *,
double *,
Modified: mc/1D/hc/trunk/ggrd_test.c
===================================================================
--- mc/1D/hc/trunk/ggrd_test.c 2008-03-31 04:21:42 UTC (rev 11650)
+++ mc/1D/hc/trunk/ggrd_test.c 2008-03-31 04:36:28 UTC (rev 11651)
@@ -7,51 +7,63 @@
static int order = 3;
hc_boolean calc_derivatives ;
double lon,lat,age;
- /*
- initialize velocity structure
- */
- ggrd = (struct ggrd_master *)calloc(1,sizeof(struct ggrd_master));
- ggrd_init_master(ggrd);
- ggrd->v.history = TRUE; /* expect history */
- ggrd->age_control = TRUE; /* expect seafloor age files */
- ggrd->age_bandlim = 900;
- /*
- read in velocities
- */
- if(ggrd_read_vel_grids(ggrd,1.0,FALSE,TRUE,
- "/home/walter/becker/data/plates/past/clb/hall/")){
- fprintf(stderr,"error opening grids\n");
- exit(-1);
- }
+ int mode;
- if(argc>1)
- sscanf(argv[1],"%lf",&time);
- else
- time = 0.0;
- dtrange = 1.0; /* transition width, in Ma */
+ mode = 1;
- fprintf(stderr,"%s: using time %g\n",argv[0],time);
- calc_derivatives = FALSE;
-
- xloc[HC_R] = HC_ND_RADIUS(0.0);
- for(lat=-89;lat<=89;lat+=2)
- for(lon=0;lon<=358;lon+=2){
- //lon=270;lat=-15;{
- xloc[HC_THETA] = LAT2THETA(lat);
- xloc[HC_PHI] = LON2PHI(lon);
- /*
- interpolate
- */
- 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],
- xloc[HC_PHI],time,ggrd, &age))
- exit(-1);
-
- fprintf(stdout,"%11g %11g\t%11g %11g %11g\t%11g\n",lon,lat,vphi[0],-vtheta[0],vr[0],age);
+ switch(mode){
+ case 1:
+ if(argc>1)
+ ggrd_grdinfo(argv[1]);
+ break;
+ case 2:
+ /*
+ initialize velocity structure
+ */
+ ggrd = (struct ggrd_master *)calloc(1,sizeof(struct ggrd_master));
+ ggrd_init_master(ggrd);
+ ggrd->v.history = TRUE; /* expect history */
+ ggrd->age_control = TRUE; /* expect seafloor age files */
+ ggrd->age_bandlim = 900;
+ /*
+ read in velocities
+ */
+ if(ggrd_read_vel_grids(ggrd,1.0,FALSE,TRUE,
+ "/home/walter/becker/data/plates/past/clb/hall/")){
+ fprintf(stderr,"error opening grids\n");
+ exit(-1);
}
-
+
+ if(argc>1)
+ sscanf(argv[1],"%lf",&time);
+ else
+ time = 0.0;
+ dtrange = 1.0; /* transition width, in Ma */
+
+ fprintf(stderr,"%s: using time %g\n",argv[0],time);
+
+ calc_derivatives = FALSE;
+
+ xloc[HC_R] = HC_ND_RADIUS(0.0);
+ for(lat=-89;lat<=89;lat+=2)
+ for(lon=0;lon<=358;lon+=2){
+ //lon=270;lat=-15;{
+ xloc[HC_THETA] = LAT2THETA(lat);
+ xloc[HC_PHI] = LON2PHI(lon);
+ /*
+ interpolate
+ */
+ 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],
+ xloc[HC_PHI],time,ggrd, &age))
+ exit(-1);
+
+ fprintf(stdout,"%11g %11g\t%11g %11g %11g\t%11g\n",lon,lat,vphi[0],-vtheta[0],vr[0],age);
+ }
+ break; /* end mode 2 */
+ }
return 0;
}
Modified: mc/1D/hc/trunk/hc.h
===================================================================
--- mc/1D/hc/trunk/hc.h 2008-03-31 04:21:42 UTC (rev 11650)
+++ mc/1D/hc/trunk/hc.h 2008-03-31 04:36:28 UTC (rev 11651)
@@ -8,14 +8,16 @@
*/
#ifndef __HC_READ_HEADER_FILE__
+#include "hc_filenames.h"
+#include <math.h>
#ifndef __GMT_INCLUDED__
+#define _GMT_MATH_H /* gmt_math.h led to clashes with some other codes */
#include "gmt.h"
-#include "gmt_bcr.h"
#define __GMT_INCLUDED__
#endif
-#include "hc_filenames.h"
+
/*
general variable type defines
@@ -25,8 +27,8 @@
#define hc_boolean unsigned short
#endif
#ifndef HC_PREC /*
- precision for most C functions
- */
+ precision for most C functions
+ */
#define HC_PREC double
#define HC_FLT_FORMAT "%lf"
#define HC_TWO_FLT_FORMAT "%lf %lf"
@@ -65,6 +67,9 @@
*/
#include "prem.h"
+
+
+
/*
dealing with velocity grids
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2008-03-31 04:21:42 UTC (rev 11650)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2008-03-31 04:36:28 UTC (rev 11651)
@@ -1,5 +1,6 @@
/* ggrd_grdtrack_util.c */
void ggrd_init_master(struct ggrd_master *);
+void ggrd_grdinfo(char *);
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);
@@ -8,7 +9,7 @@
unsigned char ggrd_grdtrack_interpolate_xy(double, double, struct ggrd_gt *, double *, unsigned char);
void ggrd_grdtrack_free_gstruc(struct ggrd_gt *);
void ggrd_find_spherical_vel_from_rigid_cart_rot(double *, double *, double *, double *, double *);
-void ggrd_print_layer_avg(float *, float *, int, int, FILE *);
+void ggrd_print_layer_avg(float *, float *, int, int, int, FILE *, int *);
int ggrd_init_thist_from_file(struct ggrd_t *, char *, unsigned char, unsigned char);
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);
@@ -128,6 +129,8 @@
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);
+/* rotvec2vel.c */
+FILE *myopen(const char *, const char *);
/* sh_ana.c */
/* shana_sh.c */
void shana_compute_allplm(int, int, double *, double *, struct shana_module *);
Modified: mc/1D/hc/trunk/make_tar
===================================================================
--- mc/1D/hc/trunk/make_tar 2008-03-31 04:21:42 UTC (rev 11650)
+++ mc/1D/hc/trunk/make_tar 2008-03-31 04:36:28 UTC (rev 11651)
@@ -2,7 +2,7 @@
#
# create a tar file
#
-ver=${1-0.1.3}
+ver=${1-0.1.4}
date=`date '+%m%d%y'`
tf=$HOME/tmp/hc-$ver.$date.tgz
More information about the cig-commits
mailing list