/* part of the shansyn spherical harmonics package, see COPYRIGHT for license */
/* $Id: spear.c,v 1.5 2001/10/10 14:29:45 becker Exp $ */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "function_macros.h"
#include "trig_constants.h"
#include "precision.h"
#include "spherical_harmonics_functions.h"
#include "legendre_macros.h"
#include "spear.h"

//
// calculate spearman rank correlation
//
COMP_PRECISION spearman_corr(COMP_PRECISION *a, 
			     COMP_PRECISION *b,
			     COMP_PRECISION *c, 
			     COMP_PRECISION *d,
			     int l)
{
  COMP_PRECISION *data1,*data2,dval,zd,probd,rs,probrs;
  int m,lmin,lmax;
  unsigned long n;
  if(l<0){// sum over all l if l given negative
    lmin=  0;
    lmax= -l;
  }else{
    lmin = l;
    lmax = l;
  }
  n=0;
  data1=calloc(sizeof(COMP_PRECISION),1);
  data2=calloc(sizeof(COMP_PRECISION),1);
  for(l=lmin;l <= lmax;l++)
    for(m=0;m<=l;m++){
      n++;
      if((data1=(COMP_PRECISION *)realloc(data1,sizeof(COMP_PRECISION)*n))==NULL)MEMERROR;
      if((data2=(COMP_PRECISION *)realloc(data2,sizeof(COMP_PRECISION)*n))==NULL)MEMERROR;
      data1[n]=   *(a+POSLM(l, m));
      data2[n]=   *(c+POSLM(l, m));
      if(m!=0){
	n++;
	if((data1=(COMP_PRECISION *)realloc(data1,sizeof(COMP_PRECISION)*n))==NULL)MEMERROR;
	if((data2=(COMP_PRECISION *)realloc(data2,sizeof(COMP_PRECISION)*n))==NULL)MEMERROR;
	data1[n+1]= *(b+POSLM(l, m));
	data2[n+1]= *(d+POSLM(l, m));
      }
    }
  spear(data1-1,data2-1,n,&dval,&zd,&probd,&rs,&probrs);
  free(data1);
  free(data2);
  return(rs);
}
/*
  
  spear routine from numerical recipes

 */
void spear(COMP_PRECISION *data1,COMP_PRECISION *data2,unsigned long n,
	   COMP_PRECISION *d,COMP_PRECISION *zd,COMP_PRECISION *probd,COMP_PRECISION *rs,COMP_PRECISION *probrs)
{
  unsigned long j;
  COMP_PRECISION vard,t,sg,sf,fac,en3n,en,df,aved,*wksp1,*wksp2,tmpd;
  
  wksp1=vector(1,n);
  wksp2=vector(1,n);
  for (j=1;j<=n;j++) {
    wksp1[j]=data1[j];
    wksp2[j]=data2[j];
  }
  sort2(n,wksp1,wksp2);
  crank(n,wksp1,&sf);
  sort2(n,wksp2,wksp1);
  crank(n,wksp2,&sg);
  *d=0.0;
  for (j=1;j<=n;j++){
    tmpd = wksp1[j]-wksp2[j];
    *d += SQUARE(tmpd);
  }
  en=n;
  en3n=en*en*en-en;
  aved=en3n/6.0-(sf+sg)/12.0;
  fac=(1.0-sf/en3n)*(1.0-sg/en3n);
  vard=((en-1.0)*en*en*SQUARE(en+1.0)/36.0)*fac;
  *zd=(*d-aved)/sqrt(vard);
  *probd=erfcc(fabs(*zd)/1.4142136);
  *rs=(1.0-(6.0/en3n)*(*d+(sf+sg)/12.0))/sqrt(fac);
  fac=(*rs+1.0)*(1.0-(*rs));
  if (fac > 0.0) {
    t=(*rs)*sqrt((en-2.0)/fac);
    df=en-2.0;
    *probrs=betai(0.5*df,0.5,df/(df+t*t));
  } else
    *probrs=0.0;
  free_vector(wksp2,1,n);
  free_vector(wksp1,1,n);
}

COMP_PRECISION betai(COMP_PRECISION a,COMP_PRECISION b,
		     COMP_PRECISION x)
{
  COMP_PRECISION bt;
  
  if (x < 0.0 || x > 1.0) nrerror("Bad x in routine betai");
  if (x == 0.0 || x == 1.0) bt=0.0;
  else
    bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
  if (x < (a+1.0)/(a+b+2.0))
    return bt*betacf(a,b,x)/a;
  else
    return 1.0-bt*betacf(b,a,1.0-x)/b;
}

void crank(unsigned long n,COMP_PRECISION *w,COMP_PRECISION *s)
{
  unsigned long j=1,ji,jt;
  float t,rank;
  
  *s=0.0;
  while (j < n) {
    if (w[j+1] != w[j]) {
      w[j]=j;
      ++j;
    } else {
      for (jt=j+1;jt<=n && w[jt]==w[j];jt++);
      rank=0.5*(j+jt-1);
      for (ji=j;ji<=(jt-1);ji++) w[ji]=rank;
      t=jt-j;
      *s += t*t*t-t;
      j=jt;
    }
  }
  if (j == n) w[n]=n;
}


#define NUMREC_MAXIT 100
#define NUMREC_EPS 3.0e-7
#define NUMREC_FPMIN 1.0e-30

COMP_PRECISION betacf(COMP_PRECISION a,COMP_PRECISION b,COMP_PRECISION x)
{
  int m,m2;
  COMP_PRECISION aa,c,d,del,h,qab,qam,qap;
  
  qab=a+b;
  qap=a+1.0;
  qam=a-1.0;
  c=1.0;
  d=1.0-qab*x/qap;
  if (fabs(d) < NUMREC_FPMIN) d=NUMREC_FPMIN;
  d=1.0/d;
  h=d;
  for (m=1;m<=NUMREC_MAXIT;m++) {
    m2=2*m;
    aa=m*(b-m)*x/((qam+m2)*(a+m2));
    d=1.0+aa*d;
    if (fabs(d) < NUMREC_FPMIN) d=NUMREC_FPMIN;
    c=1.0+aa/c;
    if (fabs(c) < NUMREC_FPMIN) c=NUMREC_FPMIN;
    d=1.0/d;
    h *= d*c;
    aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
    d=1.0+aa*d;
    if (fabs(d) < NUMREC_FPMIN) d=NUMREC_FPMIN;
    c=1.0+aa/c;
    if (fabs(c) < NUMREC_FPMIN) c=NUMREC_FPMIN;
    d=1.0/d;
    del=d*c;
    h *= del;
    if (fabs(del-1.0) < NUMREC_EPS) break;
  }
  if (m > NUMREC_MAXIT) nrerror("a or b too big, or NUMREC_MAXIT too small in betacf");
  return h;
}

#undef NUMREC_MAXIT
#undef NUMREC_EPS
#undef NUMREC_FPMIN

#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;
#define NUMREC_M 7
#define NUMREC_NSTACK 50

void sort2(unsigned long n,COMP_PRECISION *arr,COMP_PRECISION *brr)
{
  unsigned long i,ir=n,j,k,l=1;
  int *istack,jstack=0;
  COMP_PRECISION a,b,temp;
  
  istack=ivector(1,NUMREC_NSTACK);
  for (;;) {
    if (ir-l < NUMREC_M) {
      for (j=l+1;j<=ir;j++) {
	a=arr[j];
	b=brr[j];
	for (i=j-1;i>=1;i--) {
	  if (arr[i] <= a) break;
	  arr[i+1]=arr[i];
	  brr[i+1]=brr[i];
	}
	arr[i+1]=a;
	brr[i+1]=b;
      }
      if (!jstack) {
	free_ivector(istack,1,NUMREC_NSTACK);
	return;
      }
      ir=istack[jstack];
      l=istack[jstack-1];
      jstack -= 2;
    } else {
      k=(l+ir) >> 1;
      SWAP(arr[k],arr[l+1])
	SWAP(brr[k],brr[l+1])
	if (arr[l+1] > arr[ir]) {
	  SWAP(arr[l+1],arr[ir])
	    SWAP(brr[l+1],brr[ir])
	    }
      if (arr[l] > arr[ir]) {
	SWAP(arr[l],arr[ir])
	  SWAP(brr[l],brr[ir])
	  }
      if (arr[l+1] > arr[l]) {
	SWAP(arr[l+1],arr[l])
	  SWAP(brr[l+1],brr[l])
	  }
      i=l+1;
      j=ir;
      a=arr[l];
      b=brr[l];
      for (;;) {
	do i++; while (arr[i] < a);
	do j--; while (arr[j] > a);
	if (j < i) break;
	SWAP(arr[i],arr[j])
	  SWAP(brr[i],brr[j])
	  }
      arr[l]=arr[j];
      arr[j]=a;
      brr[l]=brr[j];
      brr[j]=b;
      jstack += 2;
      if (jstack > NUMREC_NSTACK) nrerror("NSTACK too small in sort2.");
      if (ir-i+1 >= j-l) {
	istack[jstack]=ir;
	istack[jstack-1]=i;
	ir=j-1;
      } else {
	istack[jstack]=j-1;
	istack[jstack-1]=l;
	l=i;
      }
    }
  }
}
#undef NUMREC_M
#undef NUMREC_NSTACK
#undef SWAP

COMP_PRECISION erfcc(COMP_PRECISION x)
{
  COMP_PRECISION t,z,ans;
  
  z=fabs(x);
  t=1.0/(1.0+0.5*z);
  ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
	  t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
			    t*(-0.82215223+t*0.17087277)))))))));
  return x >= 0.0 ? ans : 2.0-ans;
}

