[cig-commits] r12003 - mc/1D/hc/trunk
becker at geodynamics.org
becker at geodynamics.org
Thu May 22 12:07:20 PDT 2008
Author: becker
Date: 2008-05-22 12:07:20 -0700 (Thu, 22 May 2008)
New Revision: 12003
Added:
mc/1D/hc/trunk/itg-hc-geoid.ab
mc/1D/hc/trunk/sh_corr.c
Modified:
mc/1D/hc/trunk/Makefile
mc/1D/hc/trunk/hc_auto_proto.h
mc/1D/hc/trunk/rick_sh_c.c
mc/1D/hc/trunk/sh_exp.c
mc/1D/hc/trunk/sh_syn.c
Log:
Added correlation function
Modified: mc/1D/hc/trunk/Makefile
===================================================================
--- mc/1D/hc/trunk/Makefile 2008-05-22 18:29:40 UTC (rev 12002)
+++ mc/1D/hc/trunk/Makefile 2008-05-22 19:07:20 UTC (rev 12003)
@@ -132,7 +132,7 @@
all: dirs libs hc hc_extract_sh_layer \
- sh_syn sh_ana sh_power rotvec2vel
+ sh_syn sh_corr sh_ana sh_power rotvec2vel
libs: dirs hc_lib $(HEAL_LIBS) $(RICK_LIB)
@@ -158,6 +158,9 @@
sh_syn: $(LIBS) $(INCS) $(ODIR)/sh_syn.o
$(CC) $(LIB_FLAGS) $(ODIR)/sh_syn.o \
-o $(BDIR)/sh_syn -lhc -lrick $(HEAL_LIBS_LINKLINE) $(GGRD_LIBS_LINKLINE) -lm $(LDFLAGS)
+sh_corr: $(LIBS) $(INCS) $(ODIR)/sh_corr.o
+ $(CC) $(LIB_FLAGS) $(ODIR)/sh_corr.o \
+ -o $(BDIR)/sh_corr -lhc -lrick $(HEAL_LIBS_LINKLINE) $(GGRD_LIBS_LINKLINE) -lm $(LDFLAGS)
sh_power: $(LIBS) $(INCS) $(ODIR)/sh_power.o
$(CC) $(LIB_FLAGS) $(ODIR)/sh_power.o \
Modified: mc/1D/hc/trunk/hc_auto_proto.h
===================================================================
--- mc/1D/hc/trunk/hc_auto_proto.h 2008-05-22 18:29:40 UTC (rev 12002)
+++ mc/1D/hc/trunk/hc_auto_proto.h 2008-05-22 19:07:20 UTC (rev 12003)
@@ -117,12 +117,12 @@
void rick_four1_nr(float *, int, int);
/* rick_sh_c.c */
void rick_compute_allplm(int, int, float *, float *, struct rick_module *);
-void rick_compute_allplm_irreg(int, int, float *, float *, struct rick_module *, float *, int);
+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_irreg(float *, float *, int, int, float *, float *, struct rick_module *, float *, int, float *, int, unsigned short);
+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_irreg(float *, float *, int, float *, float *, int, float *, float *, struct rick_module *, float *, int, float *, int, unsigned short);
+void rick_shc2d_pre_reg(float *, float *, int, float *, float *, int, float *, float *, struct rick_module *, float *, int, float *, int, unsigned short);
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);
@@ -142,6 +142,7 @@
void shana_init(int, int, int *, int *, int *, struct shana_module *);
void shana_free_module(struct shana_module *, int);
void shana_plmbar1(double *, double *, int, int, double, struct shana_module *);
+/* sh_corr.c */
/* sh_exp.c */
void sh_allocate_and_init(struct sh_lms **, int, int, int, int, unsigned short, unsigned short);
void sh_init_expansion(struct sh_lms *, int, int, int, unsigned short, unsigned short);
@@ -149,6 +150,8 @@
void sh_clear_alm(struct sh_lms *);
double sh_total_power(struct sh_lms *);
void sh_compute_power_per_degree(struct sh_lms *, float *);
+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);
unsigned short sh_read_parameters_from_file(int *, int *, int *, int *, int *, double *, int *, FILE *, unsigned short, unsigned short, unsigned short);
void sh_print_coefficients_to_file(struct sh_lms *, int, FILE *, double *, unsigned short, unsigned short);
@@ -160,13 +163,13 @@
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_irreg(struct sh_lms *, int, unsigned short, float **, float *, int, float *, int, float *, unsigned short, 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_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_irreg_spatial_data_to_file(struct sh_lms *, int, float *, unsigned short, float, float *, int, float *, int, FILE *);
+void sh_print_reg_spatial_data_to_file(struct sh_lms *, int, float *, unsigned short, float, float *, int, float *, int, FILE *);
void sh_compute_plm(struct sh_lms *, int, float **, unsigned short);
-void sh_compute_plm_irreg(struct sh_lms *, int, float **, unsigned short, float *, int);
+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 *);
void sh_write_coeff(struct sh_lms *, int, int, int, unsigned short, double *);
void sh_aexp_equals_bexp_coeff(struct sh_lms *, struct sh_lms *);
Added: mc/1D/hc/trunk/itg-hc-geoid.ab
===================================================================
--- mc/1D/hc/trunk/itg-hc-geoid.ab (rev 0)
+++ mc/1D/hc/trunk/itg-hc-geoid.ab 2008-05-22 19:07:20 UTC (rev 12003)
@@ -0,0 +1,529 @@
+31 0 0.00000000e+00 1 1 0
+ 0 0
+ 0 0
+ 0 0
+ -1.0023972e+02 0.0000000e+00
+ 5.9891817e-03 -3.3284724e-02
+ 5.5032242e+01 -3.1590032e+01
+ 2.1593456e+01 0.0000000e+00
+ -4.5807012e+01 -5.5993746e+00
+ 2.0411921e+01 -1.3964699e+01
+ -1.6272945e+01 -3.1907573e+01
+ -1.0318153e+01 0.0000000e+00
+ 1.2095655e+01 1.0683629e+01
+ 7.9072799e+00 1.4945480e+01
+ -2.2353624e+01 4.5335619e+00
+ -4.2529837e+00 6.9665827e+00
+ 1.5491948e+00 0.0000000e+00
+ 1.4194953e+00 2.1289721e+00
+ 1.4710811e+01 -7.2948141e+00
+ 1.0193625e+01 4.8493715e+00
+ -6.6625854e+00 1.1236422e+00
+ -3.9437355e+00 1.5101141e+01
+ -3.3829442e+00 0.0000000e+00
+ 1.7127695e+00 -5.9811384e-01
+ 1.0975139e+00 -8.4326473e+00
+ -1.2914459e+00 -2.0195620e-01
+ -1.9406880e+00 -1.0635309e+01
+ 6.0272458e+00 1.2103227e+01
+ 2.1365764e-01 -5.3553200e+00
+ 2.0419429e+00 0.0000000e+00
+ -6.3367927e+00 -2.1460307e+00
+ 7.4539687e+00 2.0980007e+00
+ -5.6503147e+00 4.8981653e+00
+ -6.2038337e+00 -2.7987436e+00
+ -3.7172245e-02 -4.0445784e-01
+ -8.0944539e+00 3.4245521e+00
+ -3.4007898e-02 -5.4384912e-01
+ 1.1161641e+00 0.0000000e+00
+ -5.2250508e-01 -1.3287200e+00
+ 1.8051152e+00 1.4727212e+00
+ 4.3708783e-01 1.9393372e+00
+ -5.5127463e+00 1.5748442e+00
+ 5.7981540e-01 -2.0124208e+00
+ -1.4881598e+00 6.9698041e+00
+ -1.5173098e+00 -1.6890268e+00
+ -2.7979396e+00 2.7196367e+00
+ 6.3208469e-01 0.0000000e+00
+ -3.2069193e+00 -4.8279216e-01
+ 4.8310736e-01 -7.1511279e-01
+ 3.6233970e+00 1.6754301e+00
+ -2.1128020e-01 4.4900216e-01
+ 3.6802850e-01 1.2191248e+00
+ 1.4164891e+00 5.0300063e+00
+ 2.6617036e+00 2.1865551e+00
+ 4.2443317e+00 -6.7801553e-02
+ 1.0728765e+00 -2.1856121e+00
+ 1.2031287e+00 0.0000000e+00
+ -1.8896684e+00 2.9574283e+00
+ -2.1203921e+00 -1.1567510e+00
+ 1.5807937e-01 3.4773796e+00
+ -1.9056687e+00 -1.7828076e+00
+ 1.1119633e+00 1.1418402e+00
+ -8.4791117e-01 -1.7995765e+00
+ -1.8639208e-01 6.8785588e-02
+ 9.1588991e-01 -2.0690544e+00
+ -2.8284828e+00 8.5600446e-01
+ 2.2658239e+00 -5.3827078e-01
+ -1.1453289e+00 0.0000000e+00
+ -3.5222233e-01 6.1190386e-01
+ 4.5375884e-01 -2.2334377e+00
+ 6.8982169e-01 3.3577086e+00
+ -8.5614554e-01 -1.4385761e+00
+ -8.4417388e-01 -1.1187634e+00
+ -3.5290056e-02 7.7320665e-01
+ -1.0500768e-01 2.0264469e+00
+ -1.4216682e-01 5.5372484e-01
+ 7.0099889e-01 -9.4905590e-01
+ -1.1786302e+00 -4.1559002e-01
+ -1.0430348e+00 1.5717731e+00
+ 8.2199635e-01 0.0000000e+00
+ 1.2088857e+00 9.7381114e-01
+ 3.2185273e-01 7.0147079e-01
+ -8.9384810e-01 -5.6540176e-01
+ -1.5279464e+00 8.6590222e-02
+ -6.9659384e-01 -1.7124444e-01
+ 7.0707436e-02 8.7938868e-01
+ 4.2980671e-01 -8.0599407e-01
+ -5.8400082e-01 3.8207998e-01
+ -9.4559205e-01 -5.6315255e-01
+ -1.3986121e-01 6.9799932e-01
+ -2.5638217e-01 1.4405756e-01
+ -5.4679272e-02 -2.5039984e-01
+ 9.4140827e-01 0.0000000e+00
+ 1.1605279e+00 -8.7286557e-01
+ 1.2478291e+00 -1.4143773e+00
+ 4.8632424e-01 -2.2038000e+00
+ -8.2372445e-02 -2.6510731e-01
+ -1.3168257e+00 -1.5165763e+00
+ -7.9060102e-01 -1.4153100e-01
+ -6.7998200e-02 1.6515394e-01
+ -2.2679966e-01 -2.2239237e-01
+ -5.5881454e-01 -1.0352182e+00
+ 9.2739192e-01 -8.3111487e-01
+ 1.0043956e+00 1.0922263e-01
+ -7.0641913e-01 1.9838634e+00
+ 1.3806751e+00 -1.5374596e+00
+ -5.1139004e-01 0.0000000e+00
+ 4.2350515e-01 -6.5108287e-01
+ -8.1032108e-01 -9.1441287e-02
+ -8.2375264e-01 -4.4429870e-01
+ 3.6137309e-02 -5.1126369e-01
+ -6.6121361e-01 3.7876745e-01
+ -4.3015925e-01 5.5420826e-02
+ -8.4892282e-01 8.8742584e-02
+ -7.8828192e-01 -3.4849513e-01
+ -7.2082834e-01 -6.4214961e-01
+ 8.7534215e-01 -2.9181685e-02
+ -3.5300707e-01 8.8074607e-01
+ 1.9092763e-01 -7.0208861e-01
+ -7.2741311e-01 -1.0185179e+00
+ -1.1700691e+00 -1.0865232e-01
+ 4.9454888e-02 0.0000000e+00
+ -2.1272680e-01 -2.3650262e-01
+ -4.6316115e-01 -6.8358100e-01
+ -1.2050653e+00 -3.9848473e-01
+ -9.0627958e-01 1.5371019e-01
+ -2.7618838e-01 -1.7192322e-01
+ 7.4114269e-01 -8.2273894e-01
+ -1.3457894e+00 -1.1446745e-01
+ -7.2389028e-01 5.0021875e-01
+ -3.0004634e-01 -8.5708340e-01
+ 2.3155261e-01 3.3141911e-01
+ 2.9432425e-02 -4.1783286e-01
+ -7.3127190e-01 3.5211581e-01
+ 6.3991226e-01 1.0322347e-01
+ 1.1728119e-01 -5.5034952e-01
+ 4.2963873e-01 1.0601846e-01
+ -1.0626554e-01 0.0000000e+00
+ -5.9073604e-01 -7.5219992e-01
+ -5.5298298e-01 6.3238694e-01
+ 7.6511252e-01 4.8143035e-01
+ 9.2166187e-01 1.0825979e+00
+ 2.7344682e-01 7.7670129e-02
+ 3.1301226e-01 -8.0306187e-01
+ 1.8187535e-01 1.9518321e-01
+ -4.7836912e-01 1.2198521e-01
+ 5.0568423e-01 8.9492010e-01
+ -2.6635245e-01 2.6028328e-01
+ -4.3116111e-01 7.2203809e-02
+ 4.4134301e-01 1.5172437e-01
+ -3.1074854e-01 -2.3630387e-02
+ -4.3642485e-01 -8.7199011e-01
+ 3.2530330e-01 7.3942321e-01
+ -8.6402739e-01 6.6747523e-02
+ 4.3286993e-01 0.0000000e+00
+ 5.7222055e-01 7.1524402e-01
+ -4.5349296e-01 1.5371763e-01
+ -1.4231667e-01 -1.1467386e-01
+ 1.4607870e-01 5.7156509e-01
+ 3.6590775e-01 -1.8115635e-01
+ -2.6467330e-01 -6.6416580e-01
+ -5.6346591e-01 9.8937207e-02
+ 8.7985021e-01 8.1991664e-02
+ -7.8486124e-02 6.2361584e-01
+ -8.5634884e-02 4.1460143e-01
+ 3.6135299e-01 -2.5012262e-01
+ 6.4796922e-01 4.6153512e-01
+ -3.7238044e-01 -4.5425784e-01
+ -3.2203733e-01 2.6111300e-01
+ -1.2493265e-01 -1.1829011e-01
+ -6.8580343e-01 8.1277571e-02
+ 7.8276203e-01 4.4836157e-01
+ 1.3758441e-01 0.0000000e+00
+ -1.6246553e-01 8.8653578e-01
+ 3.3219747e-01 2.4440002e-01
+ 1.1380511e-01 1.3195802e-01
+ 1.2323413e+00 -1.7748956e-02
+ -1.3482940e-01 -5.8971213e-01
+ 3.0630843e-01 -2.9879688e-01
+ -1.5325979e-01 -1.6847513e-01
+ 6.8805975e-01 9.8009025e-02
+ 4.4148981e-01 -8.1510375e-01
+ 1.1765627e-01 -9.5669383e-02
+ 1.5535666e-01 -4.7810298e-02
+ -6.7113371e-01 -3.7361286e-01
+ 1.4108442e-01 7.8822363e-01
+ -1.8713371e-01 -2.8951266e-01
+ 9.1297586e-01 4.5752924e-01
+ 2.2932837e-01 1.4673458e-01
+ -7.8581491e-02 -9.8731282e-02
+ 6.7469730e-02 -2.4497702e-01
+ -7.4518371e-02 0.0000000e+00
+ 2.0239779e-01 -2.6925604e-02
+ 8.0626376e-01 -5.3470531e-02
+ 1.7056005e-01 -2.4881329e-02
+ 3.5663872e-01 -1.8370028e-01
+ -2.3437283e-01 -6.1928083e-01
+ -1.0827940e-01 4.2314167e-01
+ -1.2726791e-01 1.9721659e-01
+ 6.7429486e-01 -2.2542201e-01
+ -7.4179326e-02 -1.6328539e-01
+ -7.6429640e-01 -1.7098240e-01
+ -3.6782967e-01 -2.3365219e-01
+ -5.5041660e-02 2.1333150e-01
+ 1.7013037e-01 6.4349257e-01
+ -1.0701688e-01 -2.8990524e-01
+ 3.9682170e-01 3.1745312e-01
+ -4.9056727e-01 -1.6012297e-01
+ -6.4982502e-01 3.4612830e-01
+ 7.8979754e-01 -2.1859958e-01
+ 6.1006751e-02 -1.1725524e-01
+ 4.8637174e-01 0.0000000e+00
+ -1.2560648e-01 -1.5857456e-01
+ 4.5778188e-01 3.8768334e-01
+ 1.0717464e-01 -8.7757777e-01
+ 9.6331499e-02 -5.1043398e-01
+ 2.2843510e-01 1.8634291e-01
+ 2.7437423e-01 -9.8451805e-02
+ 4.7886297e-01 1.5912056e-02
+ 1.1494619e-01 4.8589158e-02
+ -3.8773739e-01 1.5820679e-01
+ -7.2858971e-01 -1.0805600e-01
+ -3.2598869e-01 4.3204954e-01
+ -1.4590613e-01 4.0912136e-01
+ -6.1842136e-01 -1.5259148e-01
+ 2.6003127e-01 -3.2445687e-01
+ 5.8183465e-01 1.9514898e-02
+ -2.8028906e-01 -7.7480278e-03
+ -1.0159805e-01 3.0979900e-01
+ 3.4644586e-01 -2.0246371e-02
+ 6.8362682e-02 -2.4641032e-01
+ 8.4257277e-02 -2.8639602e-01
+ 1.4099411e-01 0.0000000e+00
+ 3.6453370e-01 -6.4729359e-01
+ -1.2671197e-01 9.4276753e-02
+ -4.3193234e-01 -5.1263705e-01
+ -1.1527476e-01 4.2829635e-01
+ -5.8795953e-02 -2.4751975e-02
+ -3.0567333e-01 3.0937361e-04
+ 1.4913082e-01 -5.6878143e-02
+ -3.8873542e-01 4.8476479e-02
+ -3.5771804e-01 -2.0480807e-01
+ -2.5717672e-01 -2.5728950e-02
+ -1.5553343e-01 8.0264664e-01
+ -7.5039715e-02 3.3332000e-01
+ 4.3683844e-01 -3.1443968e-01
+ 4.5717479e-01 1.6590296e-01
+ -3.9644712e-01 -2.3392280e-01
+ 1.6903425e-01 -1.5010169e-01
+ 1.6845145e-01 1.6643616e-01
+ 5.9717756e-01 -2.4978579e-01
+ 6.1368915e-01 -3.7348322e-01
+ -6.0798013e-01 3.5912864e-01
+ -1.9083341e-01 8.2414149e-02
+ -2.4313357e-01 0.0000000e+00
+ -3.5376183e-01 8.6643561e-02
+ -5.9759507e-01 -2.6274564e-02
+ -2.5591680e-01 -2.2833821e-01
+ -8.1554047e-02 4.2248959e-01
+ -2.1840226e-03 7.1852834e-03
+ 2.2931683e-01 -1.2425161e-01
+ -3.9399795e-01 -1.0710242e-01
+ -5.3013037e-01 8.6055983e-02
+ -1.5076478e-01 -1.8897630e-01
+ 1.2695265e-01 5.0424418e-01
+ 1.0863254e-01 3.9623281e-01
+ 5.4559361e-02 -1.8517891e-01
+ 3.8710795e-01 -4.4020532e-01
+ 2.4523594e-01 1.8784929e-01
+ -5.8128768e-01 -1.0628943e-01
+ 2.5451964e-03 -1.6443931e-01
+ -1.9326983e-01 3.2937547e-01
+ 2.3003930e-01 -3.6558349e-01
+ -3.1712794e-01 7.5853495e-02
+ -3.7839421e-01 4.4178119e-01
+ 5.7112537e-01 -5.4074045e-01
+ -2.2738773e-01 5.2678635e-02
+ -5.0216558e-01 0.0000000e+00
+ -2.0697888e-01 -3.6536242e-01
+ -3.2471852e-01 -1.0345367e-01
+ 5.4317852e-01 3.9470194e-01
+ -5.3846851e-01 1.8707528e-01
+ 1.9925269e-02 -3.7148709e-03
+ -2.6908585e-01 3.6390114e-01
+ 1.5537910e-01 4.3179082e-02
+ 1.6605406e-01 4.1138950e-03
+ -4.1018077e-02 2.8989377e-01
+ 3.8039655e-01 -5.1203021e-02
+ -2.0825166e-01 -3.2246507e-01
+ 3.7147411e-01 -2.7091279e-01
+ 2.6195827e-01 1.1742585e-01
+ 1.5694840e-01 -3.3020584e-02
+ -4.2741988e-01 8.1016011e-02
+ 1.3328128e-01 2.4988992e-01
+ 1.2325442e-01 2.9178246e-01
+ 1.9274155e-01 -3.3427379e-01
+ 1.2352707e-01 -2.4041523e-01
+ 1.7981832e-01 -1.2304166e-01
+ -3.5439809e-01 -2.6488284e-01
+ -4.0549333e-01 1.0679771e-01
+ -6.8929364e-02 2.7199446e-01
+ -1.6452126e-03 0.0000000e+00
+ 6.1598456e-02 3.4576729e-02
+ 2.1656760e-02 3.4188129e-01
+ 1.0848300e-01 2.1962796e-01
+ 1.3472063e-01 1.1480998e-01
+ 1.6057970e-01 4.8337315e-01
+ 8.1713394e-02 3.0691622e-02
+ 1.3679490e-01 -1.1092091e-01
+ 3.4758564e-01 -8.1398747e-02
+ 2.6113386e-01 4.0236800e-01
+ 2.4752535e-01 4.7955371e-01
+ -3.4256127e-01 -4.0165714e-01
+ 2.5912672e-01 -1.3938104e-01
+ 6.8177105e-02 -6.6208180e-02
+ -4.5022666e-01 -4.6167424e-02
+ -1.4336619e-01 3.5998133e-01
+ 1.9060931e-01 6.1368956e-02
+ 2.6942531e-01 1.4538231e-01
+ -8.3333538e-03 -2.3248882e-01
+ 1.0030175e-01 1.9571124e-01
+ -1.2761228e-01 1.9496716e-01
+ -1.2496063e-01 -3.1853664e-01
+ 9.0145493e-02 -8.8796495e-02
+ 1.3966278e-01 1.9898954e-01
+ 2.8557416e-01 -8.3184789e-02
+ 7.2326660e-02 0.0000000e+00
+ -1.4373774e-01 2.0657954e-01
+ 5.1250463e-01 2.1086580e-01
+ 2.1840680e-01 3.2275387e-01
+ 2.2926913e-01 6.6980423e-03
+ 2.4908786e-01 9.3074486e-03
+ 3.6714461e-01 3.6938152e-03
+ -2.1641154e-01 1.4940363e-01
+ 4.6986377e-02 4.1417322e-03
+ 6.8520951e-01 -5.1782089e-01
+ 2.0293082e-01 -9.9635726e-02
+ -4.2629007e-02 -2.3281622e-01
+ -1.7507545e-01 2.6104526e-01
+ -1.7966255e-01 2.6284573e-01
+ -4.4944362e-01 1.4802332e-01
+ 1.0029836e-01 1.6720794e-01
+ 2.9453762e-02 -2.9106668e-01
+ 3.4993429e-01 8.3573540e-02
+ 3.4685428e-02 -3.3577240e-01
+ -1.7517280e-01 -2.2123461e-01
+ -1.8013856e-01 -2.1728656e-02
+ -2.4124074e-01 -1.7505974e-01
+ -3.1897540e-01 8.7111465e-02
+ -1.9368699e-01 2.8513090e-01
+ 9.5688775e-02 -1.8943624e-01
+ -2.3673224e-01 -1.1162147e-01
+ 1.3293139e-01 0.0000000e+00
+ 3.2332651e-03 1.4999758e-01
+ -3.5781366e-02 2.5952399e-01
+ -3.4076948e-01 -9.1259290e-02
+ 4.3037046e-01 -4.5631744e-01
+ -2.9677538e-01 -1.7647508e-01
+ 2.1665360e-01 -2.3722731e-01
+ 3.5808119e-02 -1.0033638e-01
+ 7.8301220e-02 4.0180721e-02
+ 2.9542376e-01 -1.6069255e-02
+ -3.3769315e-01 -1.2771729e-01
+ 1.0439430e-01 -3.7797781e-02
+ -3.8379533e-01 5.2594330e-02
+ 3.9594115e-03 -3.4665714e-02
+ 1.7298363e-01 1.8193922e-01
+ 3.0137354e-01 -1.8719412e-01
+ 3.5772698e-02 -1.4437099e-01
+ 2.7079052e-01 -1.7859097e-01
+ -2.9186725e-01 1.1047411e-01
+ 5.1081924e-02 -8.5028631e-02
+ 1.4058175e-01 -2.7024700e-01
+ 2.0667657e-01 -4.3804324e-02
+ 2.3413166e-01 1.5938631e-01
+ -2.3097382e-02 -2.4193339e-01
+ 2.0371618e-01 3.3349291e-01
+ -8.5311078e-02 1.3101471e-02
+ 1.4846744e-02 4.6851490e-02
+ 8.1522432e-02 0.0000000e+00
+ -6.4552348e-02 -3.5808322e-02
+ 1.4891017e-02 -1.3905154e-02
+ -7.9655215e-02 -2.3648980e-01
+ -3.6079116e-02 1.8179148e-01
+ -3.2957984e-01 -3.3818447e-01
+ 1.1208214e-01 1.4465984e-01
+ 2.5113196e-01 1.5578122e-01
+ -1.8132246e-01 -1.9768678e-01
+ -6.9175865e-02 -2.5424446e-01
+ -3.0485652e-01 2.8604576e-02
+ -5.0460004e-02 2.1295484e-01
+ -2.4735211e-01 4.6697266e-02
+ 1.0313841e-01 1.0543109e-01
+ 3.4601556e-01 2.6670962e-01
+ 4.8102755e-02 -3.3188471e-02
+ 5.8026149e-02 7.1513738e-02
+ -9.9821077e-02 -1.1387424e-02
+ -4.6736487e-02 1.7771690e-01
+ 6.1562968e-03 5.4146072e-02
+ -2.8682434e-02 8.2577183e-02
+ -1.1106656e-01 1.5008048e-01
+ -1.3132786e-01 6.5216976e-02
+ 1.2004280e-01 2.6168914e-01
+ 8.6750924e-03 -5.0536236e-02
+ -2.7286580e-01 -1.3557951e-01
+ -1.5396642e-01 -5.3839884e-02
+ -1.8191150e-01 -2.4222699e-02
+ -2.2046753e-01 0.0000000e+00
+ 1.1906588e-01 -1.9669382e-01
+ -3.3860217e-01 -1.8624606e-01
+ -6.5490181e-02 -2.2974317e-01
+ -3.1881387e-02 2.5458028e-01
+ -2.7012880e-01 1.0653731e-01
+ -1.0402942e-01 2.9162905e-01
+ 5.8373230e-02 -1.2751144e-01
+ -6.9242677e-02 -3.9633788e-02
+ -2.0939708e-01 2.8641909e-01
+ -2.3975338e-01 1.7974633e-01
+ 7.9958474e-02 4.1479478e-02
+ -1.9567734e-02 2.3325886e-01
+ -4.7966117e-02 -1.3311639e-01
+ -1.8322918e-01 -2.9017098e-01
+ 2.7637078e-01 4.9314586e-02
+ -1.0465953e-01 -2.9959557e-01
+ -3.1123194e-01 1.1242824e-01
+ 1.2871703e-01 -7.5259590e-02
+ -1.4151726e-01 -5.4385851e-01
+ -3.4318254e-02 1.4661188e-01
+ -1.5290099e-01 -1.4925315e-01
+ -5.0293977e-02 -1.6491870e-01
+ -1.5541011e-01 -6.0189211e-02
+ 2.4953163e-01 -3.1321351e-01
+ -1.6616211e-01 3.9370227e-01
+ 2.7227799e-01 8.6335820e-02
+ 1.8122679e-01 -2.7241928e-02
+ 1.5372504e-01 1.5075162e-01
+ -1.1905902e-01 0.0000000e+00
+ -9.0877034e-02 3.0522796e-01
+ -5.3702906e-02 -5.5087654e-02
+ -8.6031208e-02 2.6117975e-01
+ -5.7278062e-01 5.4140053e-02
+ 1.8009665e-01 -1.3437604e-01
+ 2.6282731e-01 2.3419916e-01
+ 8.0985969e-02 7.1066207e-02
+ -4.4169236e-01 2.1609992e-01
+ -1.7402571e-03 -6.0756514e-02
+ 3.1040200e-01 1.2344204e-03
+ 1.4051817e-01 -1.5492372e-01
+ -1.7596199e-02 -1.0956378e-01
+ 3.2669770e-02 5.9492696e-02
+ -6.9442552e-02 -1.0654105e-01
+ 2.1396529e-01 1.5905426e-01
+ 3.4814164e-03 -3.3123673e-01
+ 5.1497909e-02 1.0629352e-01
+ -8.0111128e-02 -6.3664864e-02
+ 1.4013413e-01 -1.5378910e-01
+ -1.8323311e-01 1.3052757e-01
+ 2.2470267e-01 1.1595139e-01
+ 2.6417959e-01 -2.7940046e-02
+ 4.3123471e-02 -5.0036191e-02
+ 7.0571105e-03 -5.0518762e-02
+ -1.3068963e-01 -1.9722821e-01
+ 1.7744805e-01 -1.5374501e-01
+ 1.7444248e-01 1.5958665e-02
+ 2.1509425e-01 -1.2944934e-01
+ -2.8886588e-01 1.1759543e-01
+ 1.4113935e-01 0.0000000e+00
+ -6.1582863e-03 -2.4524286e-02
+ -2.2474494e-01 -2.6472235e-02
+ -8.8403096e-02 3.0988028e-01
+ -3.7579355e-02 -8.4774841e-02
+ 5.9392753e-02 1.2800553e-01
+ -3.5905833e-03 1.8871242e-02
+ -5.3980003e-02 -8.6017562e-02
+ 1.0502578e-01 7.3781576e-02
+ 1.8096899e-01 2.6415468e-01
+ 1.1029435e-01 -1.1181426e-01
+ 2.1544295e-01 -2.3435416e-01
+ 3.8573832e-01 -2.4220227e-01
+ -3.1118319e-01 -6.7210073e-02
+ 1.0908323e-01 1.8323484e-01
+ -1.4788133e-02 2.3984980e-02
+ -2.3836302e-01 9.3467198e-02
+ 1.5533124e-01 1.3259762e-01
+ -2.5182574e-01 -1.7976165e-01
+ 3.0754650e-01 -6.2032192e-02
+ -1.1781631e-01 3.0094699e-01
+ 2.6308474e-01 1.2508058e-01
+ -1.1828106e-01 -1.9646229e-01
+ -1.2165449e-01 2.4759372e-01
+ -6.1759899e-02 -5.8156754e-02
+ -6.9418570e-02 3.5465357e-01
+ 3.0582871e-02 2.8314483e-01
+ 1.7437407e-01 -2.8319778e-01
+ -1.3140580e-01 -1.8257730e-01
+ -8.8850098e-02 -4.2339783e-02
+ 5.8401058e-02 1.9117275e-01
+ 1.5249474e-01 0.0000000e+00
+ -1.6064579e-01 4.1683964e-01
+ 1.4125958e-01 1.1543921e-01
+ 1.4879338e-01 1.7017073e-01
+ 2.5697974e-01 -1.1677588e-01
+ 2.5902136e-01 -1.2310657e-01
+ -1.9199188e-03 -6.6651293e-03
+ 3.3998968e-02 3.0275241e-02
+ -7.0679608e-02 4.8273408e-02
+ -2.0426367e-02 -1.3195491e-01
+ 9.8577625e-02 -1.1037145e-01
+ -4.4809704e-02 -4.8964954e-01
+ 1.2937246e-02 9.0143222e-02
+ -2.2011509e-01 -1.0528021e-01
+ -1.8172858e-01 7.2961576e-02
+ -9.3799749e-02 5.3325374e-02
+ -1.5519934e-01 1.3860747e-01
+ 6.6421023e-02 -1.5462027e-01
+ -1.1933248e-02 -5.4750100e-02
+ -4.8277693e-02 -1.0759573e-01
+ -5.5250452e-02 1.2934224e-01
+ 2.2460069e-01 -1.6764788e-01
+ -2.2312880e-01 -2.5331810e-01
+ -1.8191027e-01 -9.7109301e-02
+ -6.5239752e-02 -8.3825634e-02
+ 3.7532804e-01 3.9764376e-02
+ -2.8987244e-01 4.0198286e-02
+ 2.5658324e-02 -2.4896017e-01
+ 2.3039758e-01 6.4204623e-02
+ 4.5649074e-02 5.2468538e-02
+ -1.8620312e-02 -1.8352891e-01
+ 1.9910140e-01 3.5638735e-02
Modified: mc/1D/hc/trunk/rick_sh_c.c
===================================================================
--- mc/1D/hc/trunk/rick_sh_c.c 2008-05-22 18:29:40 UTC (rev 12002)
+++ mc/1D/hc/trunk/rick_sh_c.c 2008-05-22 19:07:20 UTC (rev 12003)
@@ -27,9 +27,9 @@
// theta array their derviatives with respect to theta, if ivec is set
// to 1 */
-void rick_compute_allplm_irreg(int lmax,int ivec,SH_RICK_PREC *plm,
- SH_RICK_PREC *dplm, struct rick_module *rick,
- float *theta, int ntheta)
+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)
{
int i,os;
@@ -126,10 +126,10 @@
}
/*
-converts on irregular basis with locations cos(theta)[], phi[] long
+converts on regular basis with locations cos(theta)[], phi[] long
*/
-void rick_shc2d_irreg(SH_RICK_PREC *cslm,SH_RICK_PREC *dslm,
+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,
@@ -142,15 +142,15 @@
exit(-1);
}
/* allocate memory */
- hc_svecalloc(&plm,ntheta*rick->lmsize,"rick_shc2d_irreg: mem 1");
+ hc_svecalloc(&plm,ntheta*rick->lmsize,"rick_shc2d_reg: mem 1");
if(ivec)
- hc_svecalloc(&dplm,ntheta*rick->lmsize,"rick_shc2d_irreg: mem 2");
+ hc_svecalloc(&dplm,ntheta*rick->lmsize,"rick_shc2d_reg: mem 2");
//
// compute the Plm first
- rick_compute_allplm_irreg(lmax,ivec,plm,dplm,rick,theta,ntheta);
+ rick_compute_allplm_reg(lmax,ivec,plm,dplm,rick,theta,ntheta);
//
// call the precomputed subroutine
- rick_shc2d_pre_irreg(cslm,dslm,lmax,plm,dplm,ivec,rdatax,rdatay,rick,theta,ntheta,
+ rick_shc2d_pre_reg(cslm,dslm,lmax,plm,dplm,ivec,rdatax,rdatay,rick,theta,ntheta,
phi,nphi,save_sincos_fac);
/* free legendre functions if not needed */
@@ -310,15 +310,15 @@
/* //
-irregular version
+regular version
// */
-void rick_shc2d_pre_irreg(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,
- my_boolean save_sincos_fac)
+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,
+ my_boolean save_sincos_fac)
{
/* //
// Legendre functions are precomputed
@@ -327,7 +327,7 @@
float *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_irreg: error: initialize modules first\n");
+ fprintf(stderr,"rick_shc2d_pre_reg: error: initialize modules first\n");
exit(-1);
}
@@ -335,12 +335,12 @@
if(ivec){
if(!rick->vector_sh_fac_init){
- fprintf(stderr,"rick_shc2d_pre_irreg: error: vector harmonics factors not initialized\n");
+ fprintf(stderr,"rick_shc2d_pre_reg: error: vector harmonics factors not initialized\n");
exit(-1);
}
}
if((lmax+1)*(lmax+2)/2 > rick->lmsize){
- fprintf(stderr,"rick_shc2d_pre_irreg: error: lmax %i out of bounds\n",lmax);
+ fprintf(stderr,"rick_shc2d_pre_reg: error: lmax %i out of bounds\n",lmax);
exit(-1);
}
@@ -349,8 +349,8 @@
*/
if((!save_sincos_fac)||(!rick->sin_cos_saved)){
- hc_svecrealloc(&rick->sfac,nphi*lmaxp1,"rick_shc2d_pre_irreg");
- hc_svecrealloc(&rick->cfac,nphi*lmaxp1,"rick_shc2d_pre_irreg");
+ hc_svecrealloc(&rick->sfac,nphi*lmaxp1,"rick_shc2d_pre_reg");
+ hc_svecrealloc(&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];
@@ -367,8 +367,8 @@
scalar
*/
- hc_svecrealloc(&loc_plma,ntheta*rick->lmsize,"rick_shc2d_pre_irreg 3");
- hc_svecrealloc(&loc_plmb,ntheta*rick->lmsize,"rick_shc2d_pre_irreg 4");
+ 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(k=k2=0;k < rick->lmsize;k++,k2+=2,ios1++){
loc_plma[ios1] = cslm[k2 ] * plm[ios1];
@@ -811,7 +811,7 @@
rick->old_lmax = lmax;
rick->old_ivec = ivec;
- /* for irregular expansions */
+ /* for regular expansions */
rick->cfac = rick->sfac = NULL;
/* end initial branch */
Added: mc/1D/hc/trunk/sh_corr.c
===================================================================
--- mc/1D/hc/trunk/sh_corr.c (rev 0)
+++ mc/1D/hc/trunk/sh_corr.c 2008-05-22 19:07:20 UTC (rev 12003)
@@ -0,0 +1,124 @@
+#include "hc.h"
+#include "gmt.h"
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+/*
+
+read in two sets of spherical harmonics coefficients (stdin) and
+compute linear correlation coefficients up to L = min(L1,L2), or
+limited by input parameter
+
+USAGE:
+
+for total correlation
+
+cat geoid.ab itg-hc-geoid.ab | sh_corr -1 0
+
+for correlation per degree
+
+cat geoid.ab itg-hc-geoid.ab | sh_corr -1 0 1
+
+
+Thorsten Becker (twb at usc.edu)
+
+
+*/
+
+int main(int argc, char **argv)
+{
+ int type,lmax[2],llim,shps,ilayer,nset,ivec,i,l;
+ /*
+ 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;
+ if(argc > 1){
+ if((strcmp(argv[1],"-h")==0)||(strcmp(argv[1],"--help")==0)||(strcmp(argv[1],"-help")==0))
+ argc = -1000;
+ else{
+ sscanf(argv[1],"%i",&llim);
+ if(llim == 0){
+ fprintf(stderr,"%s: error: llim should not be zero\n",argv[0]);
+ exit(-1);
+ }
+ }
+ }
+ if(argc > 2){
+ sscanf(argv[2],"%i",&i);
+ if(i)
+ short_format = TRUE;
+ }
+ if(argc > 3){
+ sscanf(argv[3],"%i",&i);
+ 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);
+ 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");
+ exit(-1);
+ }
+ if(verbose)
+ fprintf(stderr,"%s: waiting to read two spherical harmonic coefficients %s from stdin (use %s -h for help)\n",
+ argv[0],(short_format)?("(short format)"):("(long format)"),argv[0]);
+ /* allocate two expansions */
+ i=0;
+ while(sh_read_parameters_from_file(&type,(lmax+i),&shps,&ilayer,&nset,
+ &zlabel,&ivec,stdin,short_format,
+ binary,verbose)){
+ if(ivec){
+ fprintf(stderr,"%s: error, vector fields not implemented yet\n",argv[0]);
+ exit(-1);
+ }
+ if(verbose)
+ fprintf(stderr,"%s: expansion %i lmax %i type %i shps %i ...",
+ argv[0],i+1,lmax[i],type,shps);
+ /* input and init */
+ if(i==0){
+ sh_allocate_and_init(&exp1,shps,lmax[i],type,ivec,verbose,0);
+ sh_read_coefficients_from_file(exp1,shps,-1,stdin,binary,fac,verbose);
+ }else{
+ sh_allocate_and_init(&exp2,shps,lmax[i],type,ivec,verbose,0);
+ sh_read_coefficients_from_file(exp2,shps,-1,stdin,binary,fac,verbose);
+ }
+ if(verbose)
+ fprintf(stderr,"ok.\n");
+ i++;
+ }
+ if(i!=2){
+ fprintf(stderr,"%s: read error, expecting two and only two expansions\n",
+ argv[0]);
+ exit(-1);
+ }
+ 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(per_degree){
+ if(verbose)
+ fprintf(stderr,"%s: computing linear correlation per degree\n",argv[0]);
+ for(l=1;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));
+ }
+
+
+ free(exp1);free(exp2);
+
+ return 0;
+}
Modified: mc/1D/hc/trunk/sh_exp.c
===================================================================
--- mc/1D/hc/trunk/sh_exp.c 2008-05-22 18:29:40 UTC (rev 12002)
+++ mc/1D/hc/trunk/sh_exp.c 2008-05-22 19:07:20 UTC (rev 12003)
@@ -263,7 +263,58 @@
power[l] /= 2.0*((HC_CPREC)l)+1.0;
} /* end l loop */
}
+/* compute total correlation up to llim */
+HC_PREC sh_correlation(struct sh_lms *exp1, struct sh_lms *exp2, int llim)
+{
+ HC_PREC corr;
+ return sh_correlation_per_degree(exp1,exp2,1,llim);
+}
+
+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_boolean need_b;
+
+ sum[0]=sum[1]=sum[2]=0.0;
+
+ if((lmax > exp1->lmax)||(lmax > exp2->lmax)||(lmax < 1)||(lmin < 1)){
+ fprintf(stderr,"sh_compute_correlation_per_degree: error: L1 %i L2 %i lmin %i lmax %i\n",
+ exp1->lmax,exp2->lmax,lmin,lmax);
+ exit(-1);
+ }
+ for(l=lmin;l <= lmax;l++){
+
+ for(m=0;m<=l;m++){
+
+ need_b = (hc_boolean) ((m == 0) ? (0) : (2));
+ sh_get_coeff(exp1,l,m,need_b,TRUE,value1); /* convert to DT normalization */
+ sh_get_coeff(exp2,l,m,need_b,TRUE,value2); /* convert to DT normalization */
+
+ atmp = value1[0];
+ ctmp = value2[0];
+ sum[0] += atmp * ctmp;
+ sum[1] += atmp * atmp;
+ sum[2] += ctmp * ctmp;
+
+ if(need_b){
+ btmp = value1[1];
+ dtmp = value2[1];
+ sum[0] += btmp* dtmp;
+ sum[1] += btmp * btmp;
+ sum[2] += dtmp * dtmp;
+ }
+ } /* end m loop */
+
+ } /* end l loop */
+ tmp = sqrt(sum[1]*sum[2]);
+ return sum[0]/tmp;
+}
+
+
+
+
/*
print one line with all parameters needed to identify a spherical
@@ -1107,13 +1158,13 @@
/*
-compute a spatial expansion on an irregular grid given in theta and
+compute a spatial expansion on an regular grid given in theta and
phi arrays of npoints length
cos(theta) and phi are cos(colatitude) and longitude in radians
*/
-void sh_compute_spatial_irreg(struct sh_lms *exp, int ivec,
+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,
@@ -1122,7 +1173,7 @@
int npoints;
npoints = nphi * ntheta;
if((!exp[0].spectral_init)||(ivec && !exp[1].spectral_init)){
- fprintf(stderr,"sh_compute_spatial_irreg: coefficients set not initialized, ivec: %i\n",
+ fprintf(stderr,"sh_compute_spatial_reg: coefficients set not initialized, ivec: %i\n",
ivec);
exit(-1);
}
@@ -1134,33 +1185,33 @@
this routine will only compute the plm once and performs
some sanity checks
*/
- sh_compute_plm_irreg(exp,ivec,plm,verbose,theta,ntheta);
+ sh_compute_plm_reg(exp,ivec,plm,verbose,theta,ntheta);
}
switch(exp->type){
#ifdef HC_USE_HEALPIX
case SH_HEALPIX:
- HC_ERROR("sh_compute_spatial_irreg","healpix: not implemented yet");
+ HC_ERROR("sh_compute_spatial_reg","healpix: not implemented yet");
#endif
case SH_RICK:
#ifdef NO_RICK_FORTRAN
if(save_plm)
- rick_shc2d_pre_irreg(exp[0].alm,exp[1].alm,exp[0].lmax,
+ rick_shc2d_pre_reg(exp[0].alm,exp[1].alm,exp[0].lmax,
*plm,(*plm+exp->n_plm),
ivec,data,(data+npoints),
&exp->rick,theta,ntheta,phi,nphi,
save_sincos_fac);
else
- rick_shc2d_irreg(exp[0].alm,exp[1].alm,exp[0].lmax,ivec,
+ rick_shc2d_reg(exp[0].alm,exp[1].alm,exp[0].lmax,ivec,
data,(data+npoints),&exp->rick,
theta,ntheta,phi,nphi,save_sincos_fac);
#else
- HC_ERROR("sh_compute_spatial_irreg","Rick fortran not implemented");
+ HC_ERROR("sh_compute_spatial_reg","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");
+ HC_ERROR("sh_compute_spatial_reg","Spherepack not implemented");
break;
#endif
default:
@@ -1299,14 +1350,14 @@
/*
-irregular version
+regular grid 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,int ntheta,
- float *phi,int nphi,
- FILE *out)
+*/
+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,
+ FILE *out)
{
int i,j,k,l,npoints;
float lon,lat;
@@ -1318,7 +1369,7 @@
#ifdef HC_USE_HEALPIX
case SH_HEALPIX:
- HC_ERROR("sh_print_irreg_spatial_data_to_file","healpix not implemented");
+ HC_ERROR("sh_print_reg_spatial_data_to_file","healpix not implemented");
break;
#endif
case SH_RICK:
@@ -1344,7 +1395,7 @@
#ifdef HC_USE_SPHEREPACK
case SH_SPHEREPACK_GAUSS:
case SH_SPHEREPACK_EVEN:
- HC_ERROR("sh_print_irreg_spatial_data_to_file","spherepack not implemented");
+ HC_ERROR("sh_print_reg_spatial_data_to_file","spherepack not implemented");
break;
#endif
default:
@@ -1447,7 +1498,7 @@
/*
compute the associated Legendre functions for all (l,m) at all
-latidutinal lcoations once and only once for all irregularly placed latitudes in theya
+latidutinal lcoations once and only once for all regularly placed latitudes in theya
input:
exp: holds the expansion parameters
@@ -1458,7 +1509,7 @@
plm: will be re-allocated, has to be passed at least as NULL
*/
-void sh_compute_plm_irreg(struct sh_lms *exp,int ivec,SH_RICK_PREC **plm,
+void sh_compute_plm_reg(struct sh_lms *exp,int ivec,SH_RICK_PREC **plm,
hc_boolean verbose, float *theta, int npoints)
{
/* */
@@ -1466,8 +1517,8 @@
if((!exp->plm_computed_irr)||(exp->tn_plm_irr != exp->old_tnplm_irr)){
if((!exp->lmax)||(!exp->n_plm)||(!exp->tn_plm)){
- fprintf(stderr,"sh_compute_plm_irreg: error, expansion not initialized?\n");
- fprintf(stderr,"sh_compute_plm_irreg: lmax: %i n_plm: %i tn_plm: %i\n",
+ fprintf(stderr,"sh_compute_plm_reg: error, expansion not initialized?\n");
+ fprintf(stderr,"sh_compute_plm_reg: lmax: %i n_plm: %i tn_plm: %i\n",
exp->lmax,exp->n_plm,exp->tn_plm);
exit(-1);
}
@@ -1475,30 +1526,30 @@
allocate
*/
exp->old_tnplm_irr = exp->tn_plm_irr;
- hc_svecrealloc(plm,exp->old_tnplm_irr,"sh_compute_plm_irreg");
+ hc_svecrealloc(plm,exp->old_tnplm_irr,"sh_compute_plm_reg");
/*
compute the Legendre polynomials
*/
switch(exp->type){
#ifdef HC_USE_HEALPIX
- HC_ERROR("compute_plm_irreg","healpix not implemented");
+ HC_ERROR("compute_plm_reg","healpix not implemented");
break;
#endif
case SH_RICK:
if(verbose)
- fprintf(stderr,"sh_compute_plm_irreg: Rick: computing all Plm for lmax %i and %i points\n",
+ fprintf(stderr,"sh_compute_plm_reg: Rick: computing all Plm for lmax %i and %i points\n",
exp->lmax,npoints);
#ifdef NO_RICK_FORTRAN
- rick_compute_allplm_irreg(exp->lmax,ivec,*plm,(*plm+exp->old_tnplm_irr),&exp->rick,
+ rick_compute_allplm_reg(exp->lmax,ivec,*plm,(*plm+exp->old_tnplm_irr),&exp->rick,
theta,npoints);
#else
- HC_ERROR("compute_plm_irreg","rick fortran not implemented");
+ HC_ERROR("compute_plm_reg","rick fortran not implemented");
#endif
break;
#ifdef HC_USE_SPHEREPACK
case SH_SPHEREPACK_GAUSS:
case SH_SPHEREPACK_EVEN:
- HC_ERROR("compute_plm_irreg","spherepack implemented");
+ HC_ERROR("compute_plm_reg","spherepack implemented");
break;
#endif
default:
@@ -1515,14 +1566,14 @@
first lmax
*/
if(exp->old_lmax_irr != exp->lmax){
- fprintf(stderr,"sh_compute_plm_irreg: error: lmax initially %i, now %i\n",
+ fprintf(stderr,"sh_compute_plm_reg: error: lmax initially %i, now %i\n",
exp->old_lmax_irr,exp->lmax);
exit(-1);
}
/* check if ivec was initialized if ever used */
if(ivec > exp->old_ivec_irr){
- fprintf(stderr,"sh_compute_plm_irreg: error: plm are to be saved but routine was initialized\n");
- fprintf(stderr,"sh_compute_plm_irreg: error: with ivec: %i and now we want vectors, ivec: %i\n",
+ fprintf(stderr,"sh_compute_plm_reg: error: plm are to be saved but routine was initialized\n");
+ fprintf(stderr,"sh_compute_plm_reg: error: with ivec: %i and now we want vectors, ivec: %i\n",
exp->old_ivec_irr,ivec);
exit(-1);
}
Modified: mc/1D/hc/trunk/sh_syn.c
===================================================================
--- mc/1D/hc/trunk/sh_syn.c 2008-05-22 18:29:40 UTC (rev 12002)
+++ mc/1D/hc/trunk/sh_syn.c 2008-05-22 19:07:20 UTC (rev 12003)
@@ -21,12 +21,13 @@
/*
switches
*/
- hc_boolean verbose = TRUE, short_format = FALSE ,short_format_ivec = FALSE ,binary = FALSE,
- regular_basis = FALSE;
+ hc_boolean verbose = TRUE, short_format = FALSE ,short_format_ivec = FALSE ,binary = FALSE;
+ int regular_basis = 0;
double w,e,s,n,dx,dy;
/* */
+ FILE *in;
float *data,*theta,*phi;
- /* spacing for irregular output */
+ /* spacing for reg_ular output */
double dphi,x,y,dtheta;
HC_PREC fac[3] = {1.,1.,1.},zlabel;
SH_RICK_PREC *dummy;
@@ -49,7 +50,10 @@
}
if(argc > 3){
sscanf(argv[3],"%lf",&w);
- regular_basis = TRUE;
+ if(w == 999)
+ regular_basis = -1;
+ else
+ regular_basis = 1;
}
if(argc > 4)
sscanf(argv[4],"%lf",&e);
@@ -71,7 +75,8 @@
fprintf(stderr,"short_ivec:\n\t0: for short format, expect AB for scalar expansion\n");
fprintf(stderr,"\t1: for short format, expect poloidal toroidal AP BP AT BT for vector expansion\n\n");
fprintf(stderr,"w,e,...\n\tif none of those are set, will use Gauss latitudes and FFT divided longitudes dependening on lmax\n");
- fprintf(stderr,"\tif set, will switch to regular spaced output with -Rw/e/s/n -Idx/dy type output\n\n");
+ fprintf(stderr,"\tif w is set to anything but 999, will switch to regular spaced output with -Rw/e/s/n -Idx/dy type output\n");
+ fprintf(stderr,"\tif w is set to 999, will read lon lat in deg from \"tmp.lonlat\", and expand on those locations\n\n");
fprintf(stderr,"The output format will depend on the type of SH input.\n");
fprintf(stderr,"\tfor scalara: lon lat scalar if a single SH is read in, else lon lat zlabel scalar.\n");
fprintf(stderr,"\tfor vectors: lon lat v_theta v_phi if a single SH is read in, else lon lat zlabel v_theta v_phi.\n\n\n");
@@ -92,11 +97,11 @@
argv[0],lmax,ivec,zlabel);
/* input and init */
- sh_allocate_and_init(&exp,shps,lmax,type,ivec,verbose,((regular_basis)?(1):(0)));
+ sh_allocate_and_init(&exp,shps,lmax,type,ivec,verbose,((regular_basis != 0)?(1):(0)));
sh_read_coefficients_from_file(exp,shps,-1,stdin,binary,fac,verbose);
- if(regular_basis){
+ if(regular_basis == 1){
/*
- irregular basis output
+ regular basis output on regular grid
*/
if(verbose)
fprintf(stderr,"sh_syn: using regular spaced grid with -R%g/%g/%g/%g -I%g/%g spacing\n",
@@ -128,17 +133,40 @@
theta[j] = y;
hc_svecalloc(&data,npoints * shps,"sh_shsyn data");
/* compute the expansion */
- sh_compute_spatial_irreg(exp,ivec,FALSE,&dummy,
+ sh_compute_spatial_reg(exp,ivec,FALSE,&dummy,
theta,ntheta,phi,nphi,data,
verbose,FALSE);
/* output */
- sh_print_irreg_spatial_data_to_file(exp,shps,data,
+ sh_print_reg_spatial_data_to_file(exp,shps,data,
(nset>1)?(TRUE):(FALSE),
zlabel, theta,ntheta,
phi,nphi,stdout);
+ }else if(regular_basis == -1){
+ /* output on locations input lon lat file */
+ 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");
+ 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){
+ phi[npoints] = LON2PHI(dphi);
+ theta[npoints] = LAT2THETA(dtheta);
+ npoints++;
+ hc_svecrealloc(&phi,npoints+1,"sh_syn");
+ hc_svecrealloc(&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);
+ fprintf(stderr,"ERROR: not implemented yet\n");
+ exit(-1);
+
-
-
}else{ /* use the built in spatial basis (Gaussian) */
/* expansion */
hc_svecalloc(&data,exp[0].npoints * shps,"sh_syn");
@@ -148,8 +176,6 @@
(nset>1)?(TRUE):(FALSE),
zlabel,stdout);
}
-
-
free(exp);free(data);
}
More information about the cig-commits
mailing list