From 05c03b5203ce91e39634a09138ace955139642a1 Mon Sep 17 00:00:00 2001 From: thwbecker Date: Sun, 20 Mar 2022 10:05:35 -0500 Subject: [PATCH] Added grdgrd2correlation to this distro. --- fitxyee.h | 85 ++++ fitxyee_util.c | 889 ++++++++++++++++++++++++++++++++++++++++ grdgrd2correlation.c | 403 ++++++++++++++++++ test_grdgrd2correlation | 18 + 4 files changed, 1395 insertions(+) create mode 100644 fitxyee.h create mode 100644 fitxyee_util.c create mode 100644 grdgrd2correlation.c create mode 100755 test_grdgrd2correlation diff --git a/fitxyee.h b/fitxyee.h new file mode 100644 index 0000000..14aad61 --- /dev/null +++ b/fitxyee.h @@ -0,0 +1,85 @@ +/* fitxyee.c */ +#include +#include +#include +#include +#define PRECISION double +//#define NR_ACC 1e-7 +#define NR_ACC 1e-15 +#define NR_EPS 5e-15 +#define FOURFORMAT "%lf %lf %lf %lf" +#define TWOFORMAT "%lf %lf" + +/* + + fit a line through x y data with uncertainties in both x and y + + $Id: fitxyee.c,v 1.2 2005/07/04 16:47:43 becker Exp becker $ + + From Numerical Recipes in C, p.668 + + + major modifications: + + - using data structure instead of x,y,sigx,sigy + - removed global variables and put those into fit structure + +*/ +static PRECISION _tmp_sqrarg; +#define SQUARE(a) ((((_tmp_sqrarg=(a))) == 0.0) ? (0.0) : (_tmp_sqrarg*_tmp_sqrarg)) +#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) +static PRECISION _tmp_maxarg1,_tmp_maxarg2; +#define FMAX(a,b) (_tmp_maxarg1=(a),_tmp_maxarg2=(b),(_tmp_maxarg1) > (_tmp_maxarg2) ?\ + (_tmp_maxarg1) : (_tmp_maxarg2)) +/* + data structure +*/ +struct dp{ + PRECISION x,y,sigx,sigy; +}; +/* + fit structure +*/ +struct fits{ + int nn; + PRECISION *xx,*yy,*sx,*sy,*ww,aa,offs; +}; + + + +void fit(struct dp *, int, PRECISION *, PRECISION *, PRECISION *, PRECISION *, PRECISION *, PRECISION *); +void fitexy(struct dp *, int, PRECISION *, PRECISION *, PRECISION *, PRECISION *, PRECISION *, PRECISION *); +PRECISION brent(PRECISION, PRECISION, PRECISION, PRECISION (*)(void), struct fits *, PRECISION, PRECISION *); +PRECISION chixy(PRECISION, struct fits *); +PRECISION gammq(PRECISION, PRECISION); +void gcf(PRECISION *, PRECISION, PRECISION, PRECISION *); +void gser(PRECISION *, PRECISION, PRECISION, PRECISION *); +void nrerror(char []); +PRECISION *vector(long, long); +int *ivector(long, long); +unsigned char *cvector(long, long); +unsigned long *lvector(long, long); +PRECISION *dvector(long, long); +PRECISION **matrix(long, long, long, long); +PRECISION **dmatrix(long, long, long, long); +int **imatrix(long, long, long, long); +PRECISION **submatrix(PRECISION **, long, long, long, long, long, long); +PRECISION **convert_matrix(PRECISION *, long, long, long, long); +PRECISION ***f3tensor(long, long, long, long, long, long); +void free_vector(PRECISION *, long, long); +void free_ivector(int *, long, long); +void free_cvector(unsigned char *, long, long); +void free_lvector(unsigned long *, long, long); +void free_dvector(PRECISION *, long, long); +void free_matrix(PRECISION **, long, long, long, long); +void free_dmatrix(PRECISION **, long, long, long, long); +void free_imatrix(int **, long, long, long, long); +void free_submatrix(PRECISION **, long, long, long, long); +void free_convert_matrix(PRECISION **, long, long, long, long); +void free_f3tensor(PRECISION ***, long, long, long, long, long, long); +void avevar(struct dp *, unsigned long, PRECISION *, PRECISION *); +void fitline(PRECISION *, PRECISION *, int, PRECISION *, int, PRECISION *, PRECISION *, PRECISION *, PRECISION *, PRECISION *, PRECISION *); +void mnbrak(PRECISION *, PRECISION *, PRECISION *, PRECISION *, PRECISION *, PRECISION *, PRECISION (*)(void), struct fits *); +PRECISION zbrent(PRECISION (*)(void), struct fits *, PRECISION, PRECISION, PRECISION); +PRECISION gammln(PRECISION); +int comparef(struct dp *, struct dp *); diff --git a/fitxyee_util.c b/fitxyee_util.c new file mode 100644 index 0000000..f7c91f3 --- /dev/null +++ b/fitxyee_util.c @@ -0,0 +1,889 @@ +#include "fitxyee.h" + + +void fit(struct dp *data,int ndata,PRECISION *a,PRECISION *b, + PRECISION *siga,PRECISION *sigb,PRECISION *chi2, + PRECISION *q) +{ + + int i; + PRECISION wt,t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss,sigdat; + + *b=0.0; + + ss=0.0; + for (i=1;i<=ndata;i++) { + wt=1.0/SQUARE(data[i].sigy); + ss += wt; + sx += data[i].x*wt; + sy += data[i].y*wt; + } + + + sxoss=sx/ss; + for (i=1;i<=ndata;i++) { + t=(data[i].x-sxoss)/data[i].sigy; + st2 += t*t; + *b += t*data[i].y/data[i].sigy; + } + + *b /= st2; + *a=(sy-sx*(*b))/ss; + *siga=sqrt((1.0+sx*sx/(ss*st2))/ss); + *sigb=sqrt(1.0/st2); + *chi2=0.0; + for (i=1;i<=ndata;i++) + *chi2 += SQUARE(data[i].y-(*a)-(*b)*data[i].x); + *q=1.0; + sigdat=sqrt((*chi2)/(ndata-2)); + *siga *= sigdat; + *sigb *= sigdat; +} + + +#define POTN 1.571000 +#define BIG 1.0e30 +#define PI 3.141592653589793238462643383 +#define ACC NR_ACC + + + + +void fitexy(struct dp *data,int ndat,PRECISION *a,PRECISION *b, + PRECISION *siga,PRECISION *sigb,PRECISION *chi2, + PRECISION *q) +{ + int j; + PRECISION swap,amx,amn,var[3],varx,vary,ang[7],ch[7],scale, + bmn,bmx,d1,d2,r2,avea[3], + dum1,dum2,dum3,dum4,dum5; + struct fits fit[1]; + fit->xx=vector(1,ndat); + fit->yy=vector(1,ndat); + fit->sx=vector(1,ndat); + fit->sy=vector(1,ndat); + fit->ww=vector(1,ndat); + /* compute variance */ + avevar(data,ndat,avea,var); + varx = var[1];vary = var[2]; + /* reassign and rescale */ + scale=sqrt(varx/vary); + fit->nn=ndat; + for (j=1;j<=ndat;j++) { + //fprintf(stderr,"%i %g %g %g %g\n",j,data[j].x,data[j].y, data[j].sigx, data[j].sigy); + fit->xx[j]=data[j].x; + fit->yy[j]=data[j].y*scale; + fit->sx[j]=data[j].sigx; + fit->sy[j]=data[j].sigy*scale; + fit->ww[j]=sqrt(SQUARE(fit->sx[j])+SQUARE(fit->sy[j])); + } + fitline(fit->xx,fit->yy,fit->nn,fit->ww,1,&dum1,b,&dum2,&dum3,&dum4,&dum5); + fit->offs=ang[1]=0.0; + ang[2]=atan(*b); + ang[4]=0.0; + ang[5]=ang[2]; + ang[6]=POTN; + for (j=4;j<=6;j++) ch[j]=chixy(ang[j],fit); + mnbrak(&ang[1],&ang[2],&ang[3],&ch[1],&ch[2],&ch[3],(PRECISION (*)(void))chixy,fit); + *chi2=brent(ang[1],ang[2],ang[3],(PRECISION (*)(void))chixy,fit,ACC,b); + *chi2=chixy(*b,fit); + *a=fit->aa; + *q=gammq(0.5*(fit->nn-2),*chi2*0.5); + for (r2=0.0,j=1;j<=fit->nn;j++) r2 += fit->ww[j]; + r2=1.0/r2; + bmx=BIG; + bmn=BIG; + fit->offs=(*chi2)+1.0; + for (j=1;j<=6;j++) { + if (ch[j] > fit->offs) { + d1=fabs(ang[j]-(*b)); + while (d1 >= PI) d1 -= PI; + d2=PI-d1; + if (ang[j] < *b) { + swap=d1; + d1=d2; + d2=swap; + } + if (d1 < bmx) bmx=d1; + if (d2 < bmn) bmn=d2; + } + } + if (bmx < BIG) { + bmx=zbrent((PRECISION (*)(void))chixy,fit,*b,*b+bmx,ACC)-(*b); + amx=fit->aa-(*a); + bmn=zbrent((PRECISION (*)(void))chixy,fit,*b,*b-bmn,ACC)-(*b); + amn=fit->aa-(*a); + *sigb=sqrt(0.5*(bmx*bmx+bmn*bmn))/(scale*SQUARE(cos(*b))); + *siga=sqrt(0.5*(amx*amx+amn*amn)+r2)/scale; + } else (*sigb)=(*siga)=BIG; + *a /= scale; + *b=tan(*b)/scale; + free_vector(fit->ww,1,ndat); + free_vector(fit->sy,1,ndat); + free_vector(fit->sx,1,ndat); + free_vector(fit->yy,1,ndat); + free_vector(fit->xx,1,ndat); +} +#undef POTN +#undef BIG +#undef PI +#undef ACC + + +#define ITMAX 10000 +#define CGOLD 0.3819660 +#define ZEPS 1.0e-14 +#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); + +PRECISION brent(ax,bx,cx,f,fit,tol,xmin) +PRECISION (*f)(),*xmin,ax,bx,cx,tol; +struct fits *fit; +{ + int iter; + PRECISION a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm; + PRECISION e=0.0; + + a=(ax < cx ? ax : cx); + b=(ax > cx ? ax : cx); + x=w=v=bx; + fw=fv=fx=(*f)(x,fit); + for (iter=1;iter<=ITMAX;iter++) { + xm=0.5*(a+b); + tol2=2.0*(tol1=tol*fabs(x)+ZEPS); + if (fabs(x-xm) <= (tol2-0.5*(b-a))) { + *xmin=x; + return fx; + } + if (fabs(e) > tol1) { + r=(x-w)*(fx-fv); + q=(x-v)*(fx-fw); + p=(x-v)*q-(x-w)*r; + q=2.0*(q-r); + if (q > 0.0) p = -p; + q=fabs(q); + etemp=e; + e=d; + if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x)) + d=CGOLD*(e=(x >= xm ? a-x : b-x)); + else { + d=p/q; + u=x+d; + if (u-a < tol2 || b-u < tol2) + d=SIGN(tol1,xm-x); + } + } else { + d=CGOLD*(e=(x >= xm ? a-x : b-x)); + } + u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d)); + fu=(*f)(u,fit); + if (fu <= fx) { + if (u >= x) a=x; else b=x; + SHFT(v,w,x,u) + SHFT(fv,fw,fx,fu) + } else { + if (u < x) a=u; else b=u; + if (fu <= fw || w == x) { + v=w; + w=u; + fv=fw; + fw=fu; + } else if (fu <= fv || v == x || v == w) { + v=u; + fv=fu; + } + } + } + nrerror("Too many iterations in brent"); + *xmin=x; + return fx; +} +#undef ITMAX +#undef CGOLD +#undef ZEPS +#undef SHFT + + +#define BIG 1.0e30 + + +PRECISION chixy(bang,fit) +PRECISION bang; +struct fits *fit; +{ + int j; + PRECISION ans,avex=0.0,avey=0.0,sumw=0.0,b; + + b=tan(bang); + for (j=1;j<=fit->nn;j++) { + fit->ww[j] = SQUARE(b*fit->sx[j])+SQUARE(fit->sy[j]); + sumw += (fit->ww[j] = (fit->ww[j] == 0.0 ? BIG : 1.0/fit->ww[j])); + avex += fit->ww[j]*fit->xx[j]; + avey += fit->ww[j]*fit->yy[j]; + } + if (sumw == 0.0) sumw = BIG; + avex /= sumw; + avey /= sumw; + fit->aa=avey-b*avex; + for (ans = -(fit->offs),j=1;j<=fit->nn;j++) + ans += fit->ww[j]*SQUARE(fit->yy[j]-fit->aa-b*fit->xx[j]); + return ans; +} +#undef BIG + +PRECISION gammq(a,x) +PRECISION a,x; +{ + void gcf(),gser(); + void nrerror(); + PRECISION gamser,gammcf,gln; + + if (x < 0.0 || a <= 0.0) nrerror("Invalid arguments in routine gammq"); + if (x < (a+1.0)) { + gser(&gamser,a,x,&gln); + return 1.0-gamser; + } else { + gcf(&gammcf,a,x,&gln); + return gammcf; + } +} +#define ITMAX 10000 +#define EPS NR_EPS +#define FPMIN 1.0e-30 + +void gcf(gammcf,a,x,gln) +PRECISION *gammcf,*gln,a,x; +{ + PRECISION gammln(); + void nrerror(); + int i; + PRECISION an,b,c,d,del,h; + + *gln=gammln(a); + b=x+1.0-a; + c=1.0/FPMIN; + d=1.0/b; + h=d; + for (i=1;i<=ITMAX;i++) { + an = -i*(i-a); + b += 2.0; + d=an*d+b; + if (fabs(d) < FPMIN) d=FPMIN; + c=b+an/c; + if (fabs(c) < FPMIN) c=FPMIN; + d=1.0/d; + del=d*c; + h *= del; + if (fabs(del-1.0) < EPS) break; + } + if (i > ITMAX) { + fprintf(stderr,"%g %g %g\n",a,x,*gln); + nrerror("a too large, ITMAX too small in gcf"); + } + *gammcf=exp(-x+a*log(x)-(*gln))*h; +} +#undef ITMAX +#undef EPS +#undef FPMIN + + +#define ITMAX 1000 +#define EPS 5.0e-15 + +void gser(gamser,a,x,gln) +PRECISION *gamser,*gln,a,x; +{ + PRECISION gammln(); + void nrerror(); + int n; + PRECISION sum,del,ap; + + *gln=gammln(a); + if (x <= 0.0) { + if (x < 0.0) nrerror("x less than 0 in routine gser"); + *gamser=0.0; + return; + } else { + ap=a; + del=sum=1.0/a; + for (n=1;n<=ITMAX;n++) { + ++ap; + del *= x/ap; + sum += del; + if (fabs(del) < fabs(sum)*EPS) { + *gamser=sum*exp(-x+a*log(x)-(*gln)); + return; + } + } + nrerror("a too large, ITMAX too small in routine gser"); + return; + } +} +#undef ITMAX +#undef EPS +/* CAUTION: This is the traditional K&R C (only) version of the Numerical + Recipes utility file nrutil.c. Do not confuse this file with the + same-named file nrutil.c that is supplied in the same subdirectory or + archive as the header file nrutil.h. *That* file contains both ANSI and + traditional K&R versions, along with #ifdef macros to select the + correct version. *This* file contains only traditional K&R. */ + +#include +#define NR_END 1 +#define FREE_ARG char* + +void nrerror(error_text) +char error_text[]; +/* Numerical Recipes standard error handler */ +{ + + + fprintf(stderr,"Numerical Recipes run-time error...\n"); + fprintf(stderr,"%s\n",error_text); + fprintf(stderr,"...now exiting to system...\n"); + exit(1); +} + +PRECISION *vector(nl,nh) +long nh,nl; +/* allocate a PRECISION vector with subscript range v[nl..nh] */ +{ + PRECISION *v; + + v=(PRECISION *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(PRECISION))); + if (!v) nrerror("allocation failure in vector()"); + return v-nl+NR_END; +} + +int *ivector(nl,nh) +long nh,nl; +/* allocate an int vector with subscript range v[nl..nh] */ +{ + int *v; + + v=(int *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(int))); + if (!v) nrerror("allocation failure in ivector()"); + return v-nl+NR_END; +} + +unsigned char *cvector(nl,nh) +long nh,nl; +/* allocate an unsigned char vector with subscript range v[nl..nh] */ +{ + unsigned char *v; + + v=(unsigned char *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(unsigned char))); + if (!v) nrerror("allocation failure in cvector()"); + return v-nl+NR_END; +} + +unsigned long *lvector(nl,nh) +long nh,nl; +/* allocate an unsigned long vector with subscript range v[nl..nh] */ +{ + unsigned long *v; + + v=(unsigned long *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(long))); + if (!v) nrerror("allocation failure in lvector()"); + return v-nl+NR_END; +} + +PRECISION *dvector(nl,nh) +long nh,nl; +/* allocate a PRECISION vector with subscript range v[nl..nh] */ +{ + PRECISION *v; + + v=(PRECISION *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(PRECISION))); + if (!v) nrerror("allocation failure in dvector()"); + return v-nl+NR_END; +} + +PRECISION **matrix(nrl,nrh,ncl,nch) +long nch,ncl,nrh,nrl; +/* allocate a PRECISION matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + PRECISION **m; + + /* allocate pointers to rows */ + m=(PRECISION **) malloc((unsigned int)((nrow+NR_END)*sizeof(PRECISION*))); + if (!m) nrerror("allocation failure 1 in matrix()"); + m += NR_END; + m -= nrl; + + /* allocate rows and set pointers to them */ + m[nrl]=(PRECISION *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(PRECISION))); + if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +PRECISION **dmatrix(nrl,nrh,ncl,nch) +long nch,ncl,nrh,nrl; +/* allocate a PRECISION matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + PRECISION **m; + + /* allocate pointers to rows */ + m=(PRECISION **) malloc((unsigned int)((nrow+NR_END)*sizeof(PRECISION*))); + if (!m) nrerror("allocation failure 1 in matrix()"); + m += NR_END; + m -= nrl; + + /* allocate rows and set pointers to them */ + m[nrl]=(PRECISION *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(PRECISION))); + if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +int **imatrix(nrl,nrh,ncl,nch) +long nch,ncl,nrh,nrl; +/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */ +{ + long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; + int **m; + + /* allocate pointers to rows */ + m=(int **) malloc((unsigned int)((nrow+NR_END)*sizeof(int*))); + if (!m) nrerror("allocation failure 1 in matrix()"); + m += NR_END; + m -= nrl; + + + /* allocate rows and set pointers to them */ + m[nrl]=(int *) malloc((unsigned int)((nrow*ncol+NR_END)*sizeof(int))); + if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); + m[nrl] += NR_END; + m[nrl] -= ncl; + + for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +PRECISION **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl) +PRECISION **a; +long newcl,newrl,oldch,oldcl,oldrh,oldrl; +/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */ +{ + long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl; + PRECISION **m; + + /* allocate array of pointers to rows */ + m=(PRECISION **) malloc((unsigned int) ((nrow+NR_END)*sizeof(PRECISION*))); + if (!m) nrerror("allocation failure in submatrix()"); + m += NR_END; + m -= newrl; + + /* set pointers to rows */ + for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol; + + /* return pointer to array of pointers to rows */ + return m; +} + +PRECISION **convert_matrix(a,nrl,nrh,ncl,nch) +PRECISION *a; +long nch,ncl,nrh,nrl; +/* allocate a PRECISION matrix m[nrl..nrh][ncl..nch] that points to the matrix +declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1 +and ncol=nch-ncl+1. The routine should be called with the address +&a[0][0] as the first argument. */ +{ + long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1; + PRECISION **m; + + /* allocate pointers to rows */ + m=(PRECISION **) malloc((unsigned int) ((nrow+NR_END)*sizeof(PRECISION*))); + if (!m) nrerror("allocation failure in convert_matrix()"); + m += NR_END; + m -= nrl; + + /* set pointers to rows */ + m[nrl]=a-ncl; + for(i=1,j=nrl+1;i *fa) { + SHFT(dum,*ax,*bx,dum) + SHFT(dum,*fb,*fa,dum) + } + *cx=(*bx)+GOLD*(*bx-*ax); + *fc=(*func)(*cx,fit); + while (*fb > *fc) { + r=(*bx-*ax)*(*fb-*fc); + q=(*bx-*cx)*(*fb-*fa); + u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/ + (2.0*SIGN(FMAX(fabs(q-r),TINY),q-r)); + ulim=(*bx)+GLIMIT*(*cx-*bx); + if ((*bx-u)*(u-*cx) > 0.0) { + fu=(*func)(u,fit); + if (fu < *fc) { + *ax=(*bx); + *bx=u; + *fa=(*fb); + *fb=fu; + return; + } else if (fu > *fb) { + *cx=u; + *fc=fu; + return; + } + u=(*cx)+GOLD*(*cx-*bx); + fu=(*func)(u,fit); + } else if ((*cx-u)*(u-ulim) > 0.0) { + fu=(*func)(u,fit); + if (fu < *fc) { + SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx)) + SHFT(*fb,*fc,fu,(*func)(u,fit)) + } + } else if ((u-ulim)*(ulim-*cx) >= 0.0) { + u=ulim; + fu=(*func)(u,fit); + } else { + u=(*cx)+GOLD*(*cx-*bx); + fu=(*func)(u,fit); + } + SHFT(*ax,*bx,*cx,u) + SHFT(*fa,*fb,*fc,fu) + } +} +#undef GOLD +#undef GLIMIT +#undef TINY +#undef SHFT +#define ITMAX 1000 +#define EPS 5.0e-15 + +PRECISION zbrent(func,fit,x1,x2,tol) +PRECISION (*func)(),tol,x1,x2; +struct fits *fit; +{ + int iter; + PRECISION a=x1,b=x2,c=x2,d,e,min1,min2; + PRECISION fa=(*func)(a,fit),fb=(*func)(b,fit), + fc,p,q,r,s,tol1,xm; + + if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) + nrerror("Root must be bracketed in zbrent"); + fc=fb; + for (iter=1;iter<=ITMAX;iter++) { + if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) { + c=a; + fc=fa; + e=d=b-a; + } + if (fabs(fc) < fabs(fb)) { + a=b; + b=c; + c=a; + fa=fb; + fb=fc; + fc=fa; + } + tol1=2.0*EPS*fabs(b)+0.5*tol; + xm=0.5*(c-b); + if (fabs(xm) <= tol1 || fb == 0.0) return b; + if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) { + s=fb/fa; + if (a == c) { + p=2.0*xm*s; + q=1.0-s; + } else { + q=fa/fc; + r=fb/fc; + p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); + q=(q-1.0)*(r-1.0)*(s-1.0); + } + if (p > 0.0) q = -q; + p=fabs(p); + min1=3.0*xm*q-fabs(tol1*q); + min2=fabs(e*q); + if (2.0*p < (min1 < min2 ? min1 : min2)) { + e=d; + d=p/q; + } else { + d=xm; + e=d; + } + } else { + d=xm; + e=d; + } + a=b; + fa=fb; + if (fabs(d) > tol1) + b += d; + else + b += SIGN(tol1,xm); + fb=(*func)(b,fit); + } + nrerror("Maximum number of iterations exceeded in zbrent"); + return 0.0; +} +#undef ITMAX +#undef EPS +PRECISION gammln(xx) +PRECISION xx; +{ + double x,y,tmp,ser; + static double cof[6]={76.18009172947146,-86.50532032941677, + 24.01409824083091,-1.231739572450155, + 0.1208650973866179e-2,-0.5395239384953e-5}; + int j; + + y=x=xx; + tmp=x+5.5; + tmp -= (x+0.5)*log(tmp); + ser=1.000000000190015; + for (j=0;j<=5;j++) ser += cof[j]/++y; + return -tmp+log(2.5066282746310005*ser/x); +} + +int comparef(struct dp *a,struct dp *b) +{ + if(a->x < b->x) + return -1; + if(a->x == b->x) + return 0; + else + return 1; +} diff --git a/grdgrd2correlation.c b/grdgrd2correlation.c new file mode 100644 index 0000000..83a2bff --- /dev/null +++ b/grdgrd2correlation.c @@ -0,0 +1,403 @@ +/* + + +read in two grid files and compute the correlation of points within a +certain radius around each location + +this assumes spherical grids by default, but see switch + +grdgrd2correlation: usage: + grdgrd2correlation file1.grd file2.grd [radius, 500 km] [remove_lin, 0] [only_print_stat, 1] [is_global, 1] [lin_reg_mode, 2] + +grdgrd2correlation: output: + lon lat correlation slope number_of_points diff_from_global + + +for regional grids, this assumes 0...360 convention + + +*/ +#include "hc_ggrd.h" +#include "hc.h" +#include "fitxyee.h" + +void calc_mean(struct dp *, int ,double *); +void calc_std(struct dp *, int ,double *,double *); + +double correlation(struct dp *, double * , int ); +void calc_azi_r_lonlat(double ,double ,double ,double , + double *,double *); +double mod(double , double ); /* */ + +int main(int argc, char **argv) +{ + struct ggrd_gt g1[1],g2[1]; + char char_dummy,gmt_flags[100]; + int n,nphi,i,nnan,ret1,ret2; + + int is_global = 1; + static int use_dist_weighting = 1; + double x,y,x1,y1,dx,dy,rad_km,xloc,yloc,dri,corr,mean_ratio, + xmin,xmax,ymin,ymax,tmpx,tmpy,dxmin,dxmax, + weight_sum,weight; + double mean[2],std[2],r,a,b,siga,sigb,chi2,q,dr, + dphi,phi,aglobal,bglobal,cloc[2]; + struct dp *data; + int lin_reg_mode = 2; /* 1: errors in x 2: errors in x and y */ + int remove = 0; /* match the layer correlation? */ + int pstat = 0; /* exit after global stats */ + FILE *out; + + //double x_dump = -10, y_dump = 40; /* for sample dumps, if wanted */ + double x_dump = -999, y_dump = -999; + char out_name[500]; + + sprintf(out_name,"tmp.%g.%g.sample.dump",x_dump,y_dump); + + data=(struct dp *)malloc(sizeof(struct dp)); + + rad_km = 500; + if(argc < 3){ + fprintf(stderr,"%s: usage:\n\t%s file1.grd file2.grd [radius, %g km] [remove_lin, %i] [only_print_stat, %i] [is_global, %i] [lin_reg_mode, %i]\n\n", + argv[0],argv[0],rad_km,remove,pstat,is_global,lin_reg_mode); + fprintf(stderr,"%s: output:\n\tlon lat correlation slope number_of_points diff_from_global\n\n",argv[0]); + exit(-1); + } + if(argc > 3 ) + sscanf(argv[3],"%lf",&rad_km); + if(argc > 4 ) + sscanf(argv[4],"%i",&remove); + if(argc > 5 ) + sscanf(argv[5],"%i",&pstat); + if(argc > 6 ) + sscanf(argv[6],"%i",&is_global); + if(argc > 7 ) + sscanf(argv[7],"%i",&lin_reg_mode); + + + + if(is_global) + sprintf(gmt_flags,"-fg"); + else + sprintf(gmt_flags,"-f0x,1y"); + + fprintf(stderr,"%s: assuming grids %s and %s are %s, computing for radius %g km, remove_lin: %i, linreg mode: %i\n", + argv[0],argv[1],argv[2],(is_global)?("global"):("regional"),rad_km,remove, + lin_reg_mode); + + if(ggrd_grdtrack_init_general(FALSE,argv[1],&char_dummy, + gmt_flags,g1,TRUE,FALSE,FALSE)){ + fprintf(stderr,"%s: error reading %s\n",argv[0],argv[1]); + exit(-1); + } + if(ggrd_grdtrack_init_general(FALSE,argv[2],&char_dummy, + gmt_flags,g2,TRUE,FALSE,FALSE)){ + fprintf(stderr,"%s: error reading %s\n",argv[0],argv[2]); + exit(-1); + } + fprintf(stderr,"%s: read data ok\n",argv[0]); + /* + + compute total correlation and best-fit slopes from roughly even + area sampling - output in on different dampling + + */ + if(is_global){ + ymin = -89.75;ymax = 89.75; + //ymin = -89.5;ymax = 89.5; + xmin = 0; xmax = 360; + dy = .25; + //dy = .5; + }else{ + ymin = MAX(g1->south,g2->south); + ymax = MIN(g1->north,g2->north); + xmin = MAX(g1->west,g2->west); + xmax = MIN(g1->east,g2->east); + //dy = .2; /* output sampling is different */ + dy = .25; + } + + + + dxmin = 1e20;dxmax = -1e20; + fprintf(stderr,"%s: sampling %g-%g-%g %g-%g-%g\n",argv[0],xmin,dy/cos(((ymin+ymax)/2)/GGRD_PIF),xmax,ymin,dy,ymax); + for(n=nnan=0,y1 = ymax+1e-7,y = ymin;y <= y1;y += dy){ + dx = dy/cos(y/GGRD_PIF); /* adjust for sphericity */ + if(dx < dxmin) + dxmin = dx; + if(dx > dxmax) + dxmax = dx; + for(x=xmin,x1=xmax-dx+1e-7;x <= x1;x += dx){ + + if((!ggrd_grdtrack_interpolate_xy(x,y,g1,&tmpx,FALSE))|| + (!ggrd_grdtrack_interpolate_xy(x,y,g2,&tmpy,FALSE))){ + fprintf(stderr,"%s: interpolation error for lon %g lat %g\n",argv[0],x,y); + exit(-1); + } + + if(finite(tmpx) && finite(tmpy)){ /* only use non NaN */ + data[n].x = tmpx; + data[n].y = tmpy; + /* uncertainties? */ + data[n].sigx = data[n].sigy = 1.0; + data=(struct dp *)realloc(data,sizeof(struct dp)*(n+2)); + n++; + }else{ + nnan++; + } + } + } + fprintf(stderr,"%s: scanned region -R%g/%g/%g/%g, -I(%g-%g)/%g, nsample: %i nnan: %i\n", + argv[0],xmin,xmax,ymin, ymax,dxmin,dxmax,dy,n,nnan); + + /* + + compute global values + + */ + calc_mean(data,n,mean); + calc_std(data,n,std,mean); + corr = correlation(data,mean,n); /* compute correlation */ + //fprintf(stderr,"%s: mean: %g %g std %g %g corr %g\n",argv[0],mean[0],mean[1],std[0],std[1],corr); + switch(lin_reg_mode){ /* fit a linear relation ymod = a + b * x */ + case 1: + /* best fit slope, only error in y */ + fit((data-1),n,&aglobal, &bglobal,&siga,&sigb,&chi2,&q); + break; + case 2: + /* best fit slope */ + fitexy((data-1),n,&aglobal,&bglobal,&siga,&sigb,&chi2,&q); + break; + default: + fprintf(stderr,"%s: linear regression mode %i undefined\n",argv[0],lin_reg_mode); + exit(-1); + break; + } + + + fprintf(stderr,"%s: first grid: mean: %11g std: %11g\n",argv[0],mean[0],std[0]); + fprintf(stderr,"%s: second grid: mean: %11g std: %11g\n",argv[0],mean[1],std[1]); + fprintf(stderr,"%s: correlation: %11g best-fit: offset: %11g slope: %11g (%s)\n", + argv[0],corr,aglobal,bglobal,(lin_reg_mode==1)?("y err"):("x and y err")); + if(pstat){ + printf("%g %.5e %.5e\n",corr,aglobal,bglobal); + exit(-1); + } + /* + + regional correlation + + */ + if(remove){ + fprintf(stderr,"%s: removing trend by correcting second grid by %g + %g * x\n",argv[0],aglobal,bglobal); + } + + dri = 10; /* spacing of circle sampling */ + dr = rad_km / dri; + + if(is_global){ /* output dampling */ + ymin = -89;ymax = 89; + xmin = 0; xmax = 359; + dy = dx = 1; + }else{ + ymin = MAX(g1->south,g2->south); + ymax = MIN(g1->north,g2->north); + xmin = MAX(g1->west,g2->west); + xmax = MIN(g1->east,g2->east); + dy = dx = .2; + } + + for(y1 = ymax+1e-7,y = ymin;y <= y1;y += dy){ + for(x=xmin,x1=xmax+1e-7;x <= x1;x += dx){ + + /* start local loop */ + n = 0; + mean_ratio = 0.0; + weight_sum = 0.0; + data=(struct dp *)realloc(data,sizeof(struct dp)); + /* set up a circular sampling region around x,y */ + for(r=dr;r <= rad_km+1e-7;r += dr){ + dphi = dr / r * GGRD_TWOPI/dri; + nphi = (int)(GGRD_TWOPI/dphi); + dphi = GGRD_TWOPI/(double)nphi; + for(phi = 0;phi < GGRD_TWOPI;phi+=dphi){ + /* + compute new location + */ + calc_azi_r_lonlat(x,y,r,phi,&xloc,&yloc); + + /* + assign to data array + */ + ret1=ggrd_grdtrack_interpolate_xy(xloc,yloc,g1,&tmpx,FALSE); + ret2=ggrd_grdtrack_interpolate_xy(xloc,yloc,g2,&tmpy,FALSE); + //fprintf(stderr,"%i %i\n",ret1,ret2); + if(finite(tmpx) && finite(tmpy)){ + /* */ + data[n].x = tmpx; + data[n].y = tmpy; + if(use_dist_weighting){ + weight = 1/r; + data[n].sigx = data[n].sigy = 1/weight; + mean_ratio += log(tmpy/tmpx) * weight; + }else{ + weight = 1; + mean_ratio += log(tmpy/tmpx) * weight; + data[n].sigx = data[n].sigy = 1/weight; + } + if(remove){ + /* remove global trend from y? */ + data[n].y -= aglobal + data[n].x * bglobal; + } + /* */ + data=(struct dp *) + realloc(data,sizeof(struct dp)*(n+2)); + //fprintf(stderr,"%11g %11g %11g %11g\n",xloc,yloc,data[n].x,data[n].y); + weight_sum += weight; + n++; + } + } /* phi loop */ + } /* dr loop */ + //fprintf(stderr,"%g %g %i - %g %g\n",x,y,n,dr,rad_km); + if(x == x_dump && y == y_dump){ + fprintf(stderr,"%s: writing samples to %s\n",argv[0],out_name); + out = fopen(out_name,"w"); + for(i=0;i 5){ + /* compute means for both */ + calc_mean(data,n,mean); + calc_std(data,n,std,mean); + /* + central values + */ + calc_azi_r_lonlat(x,y,0,phi,&xloc,&yloc); + ret1=ggrd_grdtrack_interpolate_xy(xloc,yloc,g1,(cloc), FALSE); + ret2=ggrd_grdtrack_interpolate_xy(xloc,yloc,g2,(cloc+1),FALSE); + /* + output is + + lon lat correlation slope number_of_points diff_from_global mean_ratio + + */ + //fprintf(stderr,"%g %g %g %g\n",aglobal, cloc[0],bglobal,cloc[1]); + if(finite(cloc[0]) && finite(cloc[1])){ + if(std[0] > 1e-5 && std[1] > 1e-5){ + /* compute local correlation */ + corr = correlation(data,mean,n); + /* best fit with errors in x and y */ + fitexy((data-1),n,&a,&b,&siga,&sigb,&chi2,&q); + fprintf(stdout,"%11g %11g %11g %11g %5i %11g %11g\n", + x,y,corr,b,n,(aglobal + cloc[0] * bglobal)-cloc[1], + mean_ratio); + }else{ + fprintf(stdout,"%11g %11g NaN NaN %5i %11g NaN\n", + x,y,n,(aglobal + cloc[0] * bglobal)-cloc[1]); + + + } + } + } + } + } + + + +} +/* linear correlation coefficient */ +double correlation(struct dp *data, double *mean, int n) +{ + int i,nuse; + double s1,s2,s3,dx,dy; + s1 = s2 = s3 = 0.0; + for(i=0;i tmp.dat +xyz2grd tmp.dat -R0/359/-89/89 -I1 -Gtmp.r.grd + +makecpt -T-5/5/.5 -Croma -I > tmp.vs.cpt +makecpt -T-1/1/.1 -Cvik > tmp.r.cpt +grd2map tmp.1.grd tmp.vs +grd2map tmp.2.grd tmp.vs +grd2map tmp.r.grd tmp.r +epsmerge -par -x 1 -y 3 tmp.1.ps tmp.2.ps tmp.r.ps > corrt.ps +psconvert -Tf corrt.ps +rm tmp.1.ps tmp.2.ps tmp.r.ps corrt.ps tmp*grd +