[cig-commits] r12007 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Thu May 22 15:38:47 PDT 2008
Author: becker
Date: 2008-05-22 15:38:47 -0700 (Thu, 22 May 2008)
New Revision: 12007
Modified:
mc/1D/hc/trunk/README.TXT
mc/1D/hc/trunk/ggrd_grdtrack_util.c
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/rick_sh_c.c
mc/1D/hc/trunk/sh_corr.c
mc/1D/hc/trunk/sh_exp.c
mc/1D/hc/trunk/sh_syn.c
Log:
Fully implemented expansion of spherical harmonics on arbitrary locations.
Modified: mc/1D/hc/trunk/README.TXT
===================================================================
--- mc/1D/hc/trunk/README.TXT 2008-05-22 21:44:22 UTC (rev 12006)
+++ mc/1D/hc/trunk/README.TXT 2008-05-22 22:38:47 UTC (rev 12007)
@@ -237,6 +237,16 @@
cat etopo5.ab | sh_syn > etopo5.127.dat
+Note that sh_syn and sh_ana are only example implementations of the
+subroutines, there's very limited actual functionality. For a more
+useful spherical harmonics package, see shansyn at
+
+http://geodynamics.usc.edu/~becker/sdata.html
+
+
+That being said, also note helper programs sh_corr and sh_power.
+
+
SEATREE
HC is a module of the Solid Earth Teaching and Research Environment
Modified: mc/1D/hc/trunk/ggrd_grdtrack_util.c
===================================================================
--- mc/1D/hc/trunk/ggrd_grdtrack_util.c 2008-05-22 21:44:22 UTC (rev 12006)
+++ mc/1D/hc/trunk/ggrd_grdtrack_util.c 2008-05-22 22:38:47 UTC (rev 12007)
@@ -431,7 +431,7 @@
float dz1,dz2;
struct GRD_HEADER ogrd;
int i,one_or_zero,nx,ny,mx,my,nn;
- char filename[BUFSIZ*2],**cdummy;
+ char filename[BUFSIZ*2],*cdummy;
static int gmt_init = FALSE;
/*
deal with edgeinfo
@@ -542,7 +542,7 @@
fprintf(stderr,"ggrd_grdtrack_init: mem alloc ok\n");
#ifdef USE_GMT4
/* init the header */
- GMT_grd_init (*grd,0,cdummy,FALSE);
+ GMT_grd_init (*grd,0,&cdummy,FALSE);
#endif
if(*nz == 1){
if(verbose >= 2)
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2008-05-22 21:44:22 UTC (rev 12006)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2008-05-22 22:38:47 UTC (rev 12007)
@@ -123,6 +123,7 @@
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_init(int, int, int *, int *, int *, struct rick_module *, unsigned short);
@@ -164,10 +165,12 @@
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_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_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_get_coeff(struct sh_lms *, int, int, int, unsigned short, double *);
Modified: mc/1D/hc/trunk/rick_sh_c.c
===================================================================
--- mc/1D/hc/trunk/rick_sh_c.c 2008-05-22 21:44:22 UTC (rev 12006)
+++ mc/1D/hc/trunk/rick_sh_c.c 2008-05-22 22:38:47 UTC (rev 12007)
@@ -146,12 +146,12 @@
if(ivec)
hc_svecalloc(&dplm,ntheta*rick->lmsize,"rick_shc2d_reg: mem 2");
//
- // compute the Plm first
+ // compute the Plm first for all theta values
rick_compute_allplm_reg(lmax,ivec,plm,dplm,rick,theta,ntheta);
//
// call the precomputed subroutine
rick_shc2d_pre_reg(cslm,dslm,lmax,plm,dplm,ivec,rdatax,rdatay,rick,theta,ntheta,
- phi,nphi,save_sincos_fac);
+ phi,nphi,save_sincos_fac);
/* free legendre functions if not needed */
free(plm);
@@ -308,16 +308,18 @@
}
-/* //
+/*
-regular version
+regularly spaced version, data are requested on a ntheta by nphi grid,
+with nphi values in each column, located at phi[], and ntheta values
+in each row at theta[]
-// */
+*/
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,
+ struct rick_module *rick, float *theta,
+ int ntheta,float *phi,int nphi,
my_boolean save_sincos_fac)
{
/* //
@@ -330,9 +332,7 @@
fprintf(stderr,"rick_shc2d_pre_reg: error: initialize modules first\n");
exit(-1);
}
-
lmaxp1 = lmax + 1;
-
if(ivec){
if(!rick->vector_sh_fac_init){
fprintf(stderr,"rick_shc2d_pre_reg: error: vector harmonics factors not initialized\n");
@@ -369,7 +369,7 @@
*/
hc_svecrealloc(&loc_plma,ntheta*rick->lmsize,"rick_shc2d_pre_reg 3");
hc_svecrealloc(&loc_plmb,ntheta*rick->lmsize,"rick_shc2d_pre_reg 4");
- for(i=ios1=0;i < ntheta;i++){
+ 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];
loc_plmb[ios1] = cslm[k2+1] * plm[ios1];
@@ -387,7 +387,8 @@
m=0;l++;
}
rdatax[idata] += loc_plma[ios2+k] * rick->cfac[ios3+m];
- rdatax[idata] += loc_plmb[ios2+k] * rick->sfac[ios3+m];
+ if(m != 0)
+ rdatax[idata] += loc_plmb[ios2+k] * rick->sfac[ios3+m];
}
}
}
@@ -453,7 +454,120 @@
hc_svecrealloc(&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)
+{
+ /* //
+ // Legendre functions are precomputed
+ // */
+ double dpdt,dpdp,mphi,sin_theta,sfac,cfac;
+ float *loc_plma=NULL,*loc_plmb=NULL,*plm=NULL,*dplm=NULL;
+ int i,j,k,k2,m,l,lmaxp1,lm1;
+ if(!rick->initialized){
+ fprintf(stderr,"rick_shc2d_pre_reg: error: initialize modules first\n");
+ exit(-1);
+ }
+ lmaxp1 = lmax + 1;
+ if((lmax+1)*(lmax+2)/2 > rick->lmsize){
+ fprintf(stderr,"rick_shc2d_pre_reg: error: lmax %i out of bounds\n",lmax);
+ exit(-1);
+ }
+ if(ivec == 0){
+ /*
+ scalar
+
+ */
+ hc_svecrealloc(&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);
+ /* m = 0 , l = 0*/
+ l=0;m=0;
+ rdatax[i] = cslm[0] * plm[0];
+ for (k=1,k2=2; k < rick->lmsize; k++,k2+=2) {
+ m++;
+ if (m > l) {
+ m=0;l++;
+ }
+ //fprintf(stderr,"%5i %5i %11g %11g\n",l,m,cslm[k2],cslm[k2+1]);
+ if(m != 0){
+ mphi = (double)m*(double)phi[i];
+ rdatax[i] += cslm[k2] * plm[k] * cos(mphi);
+ rdatax[i] += cslm[k2+1] * plm[k] * sin(mphi);
+ }else{
+ rdatax[i] += cslm[k2] * plm[k] ;
+ }
+ }
+ } /* end data loop */
+
+ free(plm);
+ /* end scalar part */
+ } else {
+ /*
+ vector harmonics
+ */
+ hc_svecrealloc(&plm,rick->lmsize,"rick_shc2d_irreg: mem 1");
+ hc_svecrealloc(&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);
+ sin_theta = sin(theta[i]);
+
+ rdatax[i] = rdatay[i] = 0.0;
+ l = 0;m = 0;
+ /* start at l = 1 */
+ for (k=1,k2=2; k < rick->lmsize; k++,k2+=2) {
+ /* loop through l,m */
+ m++;
+ if (m > l) {
+ m=0;l++;
+ }
+ lm1 = l - 1;
+ // fprintf(stderr,"%5i %5i\t %11g %11g\tPA: %11g PB: %11g\tTA: %11g TB:%11g\n",
+ // l,m,rick->ell_factor[lm1],plm[ios2+k],
+ // SH_RICK_FACTOR(l, m) * cslm[k2] ,SH_RICK_FACTOR(l, m) *cslm[k2+1],
+ // SH_RICK_FACTOR(l, m) * dslm[k2], SH_RICK_FACTOR(l, m) *dslm[k2+1]);
+ dpdt = dplm[k] * rick->ell_factor[lm1]; /* d_theta(P_lm) factor */
+ dpdp = ((SH_RICK_PREC)m) * plm[k]/ sin_theta;
+ 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);
+ }else{
+ mphi = sfac = 0.0;
+ cfac = 1.0;
+ }
+ /*
+
+ add up contributions from all l,m
+
+ u_theta
+
+ */
+ rdatax[i] += /* cos term */
+ (cslm[k2] * dpdt + dslm[k2+1] * dpdp) * cfac;
+ rdatax[i] += /* sin term */
+ (cslm[k2+1] * dpdt - dslm[k2] * dpdp) * sfac;
+ /*
+ u_phi
+ */
+ rdatay[i] += // cos term
+ (cslm[k2+1] * dpdp - dslm[k2] * dpdt) * cfac;
+ rdatay[i] += // sin term
+ (- cslm[k2] * dpdp - dslm[k2+1] * dpdt) * sfac;
+ }
+ } /* end phi loop */
+
+ free(plm);free(dplm);
+ } /* end vector part */
+}
+
void rick_shd2c(SH_RICK_PREC *rdatax,SH_RICK_PREC *rdatay,
int lmax,int ivec,SH_RICK_PREC *cslm,
SH_RICK_PREC *dslm,struct rick_module *rick)
Modified: mc/1D/hc/trunk/sh_corr.c
===================================================================
--- mc/1D/hc/trunk/sh_corr.c 2008-05-22 21:44:22 UTC (rev 12006)
+++ mc/1D/hc/trunk/sh_corr.c 2008-05-22 22:38:47 UTC (rev 12007)
@@ -20,7 +20,11 @@
cat geoid.ab itg-hc-geoid.ab | sh_corr -1 0 1
+for correlation between l = 4 and 9
+cat geoid.ab itg-hc-geoid.ab | sh_corr 9 0 0 4
+
+
Thorsten Becker (twb at usc.edu)
@@ -28,14 +32,15 @@
int main(int argc, char **argv)
{
- int type,lmax[2],llim,shps,ilayer,nset,ivec,i,l;
+ int type,lmax[2],llim,shps,ilayer,nset,ivec,i,l,lmin;
/*
switches
*/
hc_boolean verbose = TRUE, short_format = FALSE , binary = FALSE, per_degree = FALSE;
struct sh_lms *exp1,*exp2;
HC_PREC zlabel,fac[3] = {1.,1.,1.};
- llim = -1;
+ llim = -1; /* max l */
+ lmin = 1; /* min l */
if(argc > 1){
if((strcmp(argv[1],"-h")==0)||(strcmp(argv[1],"--help")==0)||(strcmp(argv[1],"-help")==0))
argc = -1000;
@@ -57,15 +62,18 @@
if(i)
per_degree = TRUE;
}
- if((argc > 4)|| (argc < 0)){
- fprintf(stderr,"usage: cat sh.1 sh.2 | %s [llim, %i] [short_format, %i] [per_degree, %i]\n",
- argv[0],llim,short_format,per_degree);
+ if(argc > 4)
+ sscanf(argv[4],"%i",&lmin);
+ if((argc > 5)|| (argc < 0)){
+ fprintf(stderr,"usage: cat sh.1 sh.2 | %s [llim, %i] [short_format, %i] [per_degree, %i] [lmin, %i]\n",
+ argv[0],llim,short_format,per_degree,lmin);
fprintf(stderr,"reads two spherical harmonic expansions from stdin and computes\n");
fprintf(stderr,"linear correlation coefficients for L=min(L1,L2) if llim == -1\n");
fprintf(stderr,"if llim > 0, will limit the maximum degree to min(llim,L1,L2)\n");
fprintf(stderr,"short_format:\n\t0: expects regular format with long header\n");
fprintf(stderr,"\t1: expects short format with only llim in header\n");
fprintf(stderr,"per_degree: if set, will print out per degree, if not total correlation\n\n");
+ fprintf(stderr,"lmin: will compute total correlation starting at lmin\n");
exit(-1);
}
if(verbose)
@@ -100,21 +108,28 @@
argv[0]);
exit(-1);
}
+ /* check bounds */
if(llim < 0){
llim = (lmax[0] < lmax[1])?(lmax[0]):(lmax[1]);
}else{
llim = (llim < lmax[0])?(llim):(lmax[0]);
llim = (llim < lmax[1])?(llim):(lmax[1]);
}
+ if(lmin>llim){
+ fprintf(stderr,"%s: error: lmin (%i) needs to be <= llim (%i)\n",
+ argv[0],lmin,llim);
+ exit(-1);
+ }
+ /* */
if(per_degree){
if(verbose)
- fprintf(stderr,"%s: computing linear correlation per degree\n",argv[0]);
- for(l=1;l<=llim;l++)
+ 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));
}else{
if(verbose)
- fprintf(stderr,"%s: computing linear correlation up to %i\n",argv[0],llim);
- fprintf(stdout,"%14.7e\n",sh_correlation(exp1,exp2,llim));
+ 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));
}
Modified: mc/1D/hc/trunk/sh_exp.c
===================================================================
--- mc/1D/hc/trunk/sh_exp.c 2008-05-22 21:44:22 UTC (rev 12006)
+++ mc/1D/hc/trunk/sh_exp.c 2008-05-22 22:38:47 UTC (rev 12007)
@@ -945,11 +945,11 @@
}else{ /* write to file */
if(!use_3d){
/* write in lon lat format */
- fprintf(out,"%g %g\n",PHI2LON(xp[HC_PHI]),
+ fprintf(out,"%14.7f %14.7f\n",PHI2LON(xp[HC_PHI]),
THETA2LAT(xp[HC_THETA]));
}else{
/* write in lon lat z format */
- fprintf(out,"%g %g %g\n",PHI2LON(xp[HC_PHI]),
+ fprintf(out,"%14.7f %14.7f %g\n",PHI2LON(xp[HC_PHI]),
THETA2LAT(xp[HC_THETA]),z);
}
}
@@ -1219,6 +1219,40 @@
break;
}
}
+void sh_compute_spatial_irreg(struct sh_lms *exp, int ivec,
+ float *theta, float *phi,int npoints,
+ float *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",
+ ivec);
+ exit(-1);
+ }
+ switch(exp->type){
+#ifdef HC_USE_HEALPIX
+ case SH_HEALPIX:
+ HC_ERROR("sh_compute_spatial_irreg","healpix: not implemented yet");
+#endif
+ case SH_RICK:
+#ifdef NO_RICK_FORTRAN
+ rick_shc2d_irreg(exp[0].alm,exp[1].alm,exp[0].lmax,ivec,
+ data,(data+npoints),&exp->rick,
+ theta,phi,npoints);
+#else
+ HC_ERROR("sh_compute_spatial_irreg","Rick fortran not implemented");
+#endif
+ break;
+#ifdef HC_USE_SPHEREPACK
+ case SH_SPHEREPACK_GAUSS:
+ case SH_SPHEREPACK_EVEN:
+ HC_ERROR("sh_compute_spatial_irreg","Spherepack not implemented");
+ break;
+#endif
+ default:
+ sh_exp_type_error("sh_compute_spatial",exp);
+ break;
+ }
+}
/*
@@ -1247,7 +1281,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,"%12.5e ",plm[j*n_plm+i]);
+ fprintf(out,"%16.7e ",plm[j*n_plm+i]);
fprintf(out,"\n");
}
break;
@@ -1256,7 +1290,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,"%12.5e ",plm[j*n_plm+i]);
+ fprintf(out,"%16.7e ",plm[j*n_plm+i]);
fprintf(out,"\n");
}
break;
@@ -1337,13 +1371,13 @@
/* print coordinates */
if(!use_3d){
/* print lon lat */
- fprintf(out,"%11g %11g\t",lon,lat);
+ fprintf(out,"%14.7f %14.7f\t",lon,lat);
}else{
/* print lon lat z[i] */
- fprintf(out,"%11g %11g %11g\t",lon,lat,z);
+ fprintf(out,"%14.7f %14.7f %14.7f\t",lon,lat,z);
}
for(k=0;k < shps;k++) /* loop through all scalars */
- fprintf(out,"%11g ",data[j+exp[0].npoints*k]);
+ fprintf(out,"%14.7e ",data[j+exp[0].npoints*k]);
fprintf(out,"\n");
} /* end points in lateral space loop */
}
@@ -1373,20 +1407,20 @@
break;
#endif
case SH_RICK:
- for(i=l=00;i<ntheta;i++){
+ for(i=l=0;i<ntheta;i++){
lat = THETA2LAT(theta[i]);
for(j=0;j<nphi;j++,l++){
lon = PHI2LON(phi[j]);
/* print coordinates */
if(!use_3d){
/* print lon lat */
- fprintf(out,"%11g %11g\t",lon,lat);
+ fprintf(out,"%14.7f %14.7f\t",lon,lat);
}else{
/* print lon lat z[i] */
- fprintf(out,"%11g %11g %11g\t",lon,lat,z);
+ fprintf(out,"%14.7f %14.7f %14.7f\t",lon,lat,z);
}
for(k=0;k<shps;k++) /* loop through all scalars */
- fprintf(out,"%11g ",data[l+npoints*k]);
+ fprintf(out,"%14.7e ",data[l+npoints*k]);
fprintf(out,"\n");
}
}
@@ -1405,6 +1439,55 @@
}
/*
+irregular, arbitrary version
+
+*/
+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,
+ FILE *out)
+{
+ int i,k;
+ float lon,lat;
+ /*
+ get coordinates
+ */
+ switch(exp->type){
+#ifdef HC_USE_HEALPIX
+ case SH_HEALPIX:
+ HC_ERROR("sh_print_irreg_spatial_data_to_file","healpix not implemented");
+ break;
+#endif
+ case SH_RICK:
+ for(i=0;i<npoints;i++){
+ lat = THETA2LAT(theta[i]);
+ lon = PHI2LON(phi[i]);
+ /* print coordinates */
+ if(!use_3d){
+ /* print lon lat */
+ fprintf(out,"%14.7f %14.7f\t",lon,lat);
+ }else{
+ /* print lon lat z[i] */
+ fprintf(out,"%14.7f %14.7f %14.7f\t",lon,lat,z);
+ }
+ for(k=0;k < shps;k++) /* loop through all scalars */
+ fprintf(out,"%14.7e ",data[i+npoints*k]);
+ fprintf(out,"\n");
+ }
+ break;
+#ifdef HC_USE_SPHEREPACK
+ case SH_SPHEREPACK_GAUSS:
+ case SH_SPHEREPACK_EVEN:
+ HC_ERROR("sh_print_irreg_spatial_data_to_file","spherepack not implemented");
+ break;
+#endif
+ default:
+ sh_exp_type_error("sh_print_spatial_data",exp);
+ break;
+ } /* end type branch */
+}
+/*
+
compute the associated Legendre functions for all (l,m) at
all latidutinal lcoations once and only once
Modified: mc/1D/hc/trunk/sh_syn.c
===================================================================
--- mc/1D/hc/trunk/sh_syn.c 2008-05-22 21:44:22 UTC (rev 12006)
+++ mc/1D/hc/trunk/sh_syn.c 2008-05-22 22:38:47 UTC (rev 12007)
@@ -134,13 +134,13 @@
hc_svecalloc(&data,npoints * shps,"sh_shsyn data");
/* compute the expansion */
sh_compute_spatial_reg(exp,ivec,FALSE,&dummy,
- theta,ntheta,phi,nphi,data,
- verbose,FALSE);
+ theta,ntheta,phi,nphi,data,
+ verbose,FALSE);
/* output */
sh_print_reg_spatial_data_to_file(exp,shps,data,
- (nset>1)?(TRUE):(FALSE),
- zlabel, theta,ntheta,
- phi,nphi,stdout);
+ (nset>1)?(TRUE):(FALSE),
+ zlabel, theta,ntheta,
+ phi,nphi,stdout);
}else if(regular_basis == -1){
/* output on locations input lon lat file */
if(verbose)
@@ -163,11 +163,11 @@
if(verbose)
fprintf(stderr,"sh_syn: read %i locations lon lat from tmp.lonlat for expansion\n",npoints);
fclose(in);
- fprintf(stderr,"ERROR: not implemented yet\n");
- exit(-1);
-
-
- }else{ /* use the built in spatial basis (Gaussian) */
+ hc_svecalloc(&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");
sh_compute_spatial(exp,ivec,FALSE,&dummy,data,verbose);
More information about the cig-commits
mailing list