/* part of the shansyn spherical harmonics package, see COPYRIGHT for license */
/* $Id: powercorr.c,v 1.9 2001/09/19 20:21:13 becker Exp becker $ */
#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"
/*

calculates power per degree and unit area,
normalized by the 2l+1 entries


*/

COMP_PRECISION degree_power(COMP_PRECISION *a,
			    COMP_PRECISION *b, int l)
{
  COMP_PRECISION tmp;
  int m,nfac,os;
  // m=0
  tmp = SQUARE(a[POSLM(l, 0)]);
  // 1<=m<=l
  for(m=1;m<=l;m++){
    os = POSLM(l, m);
    tmp += SQUARE(a[os]);
    tmp += SQUARE(b[os]);
  }
  nfac = 2*l + 1;
  tmp /= (COMP_PRECISION)nfac;
  return tmp;
}

/*
  
  calculates power per degree and unit area for real/imaginary coeffs
  normalized by 4l+2 coefficients
  
*/
COMP_PRECISION degree_power_gsh(COMP_PRECISION *ar,
				COMP_PRECISION *ai,
				COMP_PRECISION *br,
				COMP_PRECISION *bi,
				int l)
{
  COMP_PRECISION tmp;
  int m,os,nfac;
  // m=0
  os = POSLM(l,0);
  tmp = SQUARE(ar[os]) + SQUARE(ai[os]);
  // 1<=m<=l
  for(m=1;m<=l;m++){
    os = POSLM(l, m);
    tmp += SQUARE(ar[os]) + SQUARE(ai[os]);
    tmp += SQUARE(br[os]) + SQUARE(bi[os]);
  }
  nfac = 4*l+2;
  tmp /= (COMP_PRECISION)nfac;
  return tmp;
}

/*

  calculate linear correlation coefficient between 
  two GSH 2phi or 4phi model expansions

  this sums up the dot product of all l >= 2 terms for 2phi (ialpha == 2), 
  and l >= 4 for 4phi (ialpha == 4)

*/
COMP_PRECISION correlation_gsh(COMP_PRECISION *ar, COMP_PRECISION *ai,
			       COMP_PRECISION *br, COMP_PRECISION *bi,
			       COMP_PRECISION *cr, COMP_PRECISION *ci, 
			       COMP_PRECISION *dr, COMP_PRECISION *di, 
			       int l, int weighted_for_degree,
			       int ialpha, int lmin)
{
  COMP_PRECISION sum[4],tmp,artmp,aitmp,brtmp,bitmp,w,
    crtmp,citmp,drtmp,ditmp;
  int m,lmax,os;
  sum[0]=sum[1]=sum[2]=0.0;
  if(l < 0){
    //
    // sum over all nonzero ls, lmax = -l, if l given negative
    //
    switch(ialpha){
    case 2:
      lmin = MAX(2,lmin);
      break;
    case 4:
      lmin = MAX(4,lmin);
      break;
    default:
      fprintf(stderr,"correlation_gsh: ialpha %i undefined\n",ialpha);
      exit(-1);
    }
    lmax= -l;
  }else{
    lmin = l;
    lmax = l;
  }
  for(l=lmin;l <= lmax;l++){
    if(weighted_for_degree)	/* weighted by the number of entries per degree (no good) */
      w = 1.0/(4.0*(COMP_PRECISION)l+2.0);
    else
      w = 1.0;	
    for(m=0;m <= l;m++){
      os = POSLM(l, m);
      artmp = ar[os] * w;aitmp = ai[os] * w;
      crtmp = cr[os] * w;citmp = ci[os] * w;
      sum[0] += artmp * crtmp + aitmp * citmp;
      //fprintf(stderr,"%20.4e %20.4e %20.4e %20.4e\n",artmp,crtmp,aitmp,citmp);
      sum[1] += SQUARE(artmp) + SQUARE(aitmp);
      sum[2] += SQUARE(crtmp) + SQUARE(citmp);
      if(m != 0){
	brtmp = br[os] * w;bitmp = bi[os] * w;
	drtmp = dr[os] * w;ditmp = di[os] * w;
	//fprintf(stderr,"%20.4e %20.4e %20.4e %20.4e\n",brtmp,drtmp,bitmp,ditmp);
 	sum[0] += brtmp * drtmp + bitmp * ditmp;
	sum[1] += SQUARE(brtmp) + SQUARE(bitmp);
	sum[2] += SQUARE(drtmp) + SQUARE(ditmp);
      }
    }
  }
  tmp = sqrt(sum[1] * sum[2]);
  return(sum[0]/tmp);
}
/*

calculate linear correlation coefficient between two model
expansions, weighted by each degree if weighted_by_degree is set
(this is no good)

*/
COMP_PRECISION correlation(COMP_PRECISION *a, COMP_PRECISION *b,
			   COMP_PRECISION *c, COMP_PRECISION *d, 
			   int l, int weighted_for_degree,int lmin)
{
  COMP_PRECISION sum[4],tmp,atmp,btmp,ctmp,dtmp,w;
  int m,lmax,os;
  sum[0]=sum[1]=sum[2]=0.0;
  if(l < 0){// sum over all l if l given negative
    lmin=  MAX(1,lmin);
    lmax= -l;
  }else{
    lmin = l;
    lmax = l;
  }
  for(l=lmin;l <= lmax;l++){
    if(weighted_for_degree)	/* weighted by the number of entries per degree */
      w = 1.0/(2.0*(COMP_PRECISION)l+1.0);
    else
      w = 1.0;			/* better */
    for(m=0;m<=l;m++){
      os = POSLM(l, m);
      atmp = a[os] * w;
      ctmp = c[os] * w;
      sum[0] += atmp * ctmp;
      sum[1] += SQUARE(atmp);
      sum[2] += SQUARE(ctmp);
      if(m != 0){
	btmp = b[os] * w;
	dtmp = d[os] * w;
	sum[0] += btmp* dtmp;
	sum[1] += SQUARE(btmp);
	sum[2] += SQUARE(dtmp);
      }
    }
  }
  tmp=sqrt(sum[1]*sum[2]);
  return(sum[0]/tmp);
}


COMP_PRECISION ccl_correlation(COMP_PRECISION *a, COMP_PRECISION *b,
			       int l,int lmax, int *dof)
{
  COMP_PRECISION sum[4];
  int m,os,os2;
  *dof = 0;
  if(l >= lmax){
    fprintf(stderr,"ccl_corr: can only compute r_l,l+1 up to %i-1=%i\n",
	    lmax,lmax-1);
    exit(-1);
  }
  sum[0]=sum[1]=sum[2]=0.0;
  for(m=0;m <= l;m++){
    os =  POSLM(l, m);		/* l,  m */
    os2 = POSLM(l+1, m);	/* l+1,m */
    sum[0] += a[os] * a[os2];
    sum[1] += SQUARE(a[os]);
    sum[2] += SQUARE(a[os2]);
    *dof += 1;
    if(m != 0){
      sum[0] += b[os] * b[os2];
      sum[1] += SQUARE(b[os]);
      sum[2] += SQUARE(b[os2]);
      *dof += 1;
    }
  }
  return(sum[0]/sqrt(sum[1]*sum[2]));
}


/*

  calculates correlation based on 1/(2l+1) coefficients
  (not good)

 */
COMP_PRECISION weighted_correlation(COMP_PRECISION *a, COMP_PRECISION *b,
				    COMP_PRECISION *c, COMP_PRECISION *d,int l,int lmin)
{
  return correlation(a,b,c,d,l,1,lmin);
}


/*
  calculates RMS normalized by surface area of sphere, without l = 0 term
*/
COMP_PRECISION calc_rms(COMP_PRECISION *a, COMP_PRECISION *b,int lmax)
{
  int l;
  COMP_PRECISION rms;
  rms=0.0;
  for(l=1;l<=lmax;l++)
    rms += ((COMP_PRECISION)(2*l+1))*degree_power(a,b,l);
  return(sqrt(rms)/TWO_SQRT_PI);
}
/*
  same for gsh
*/
COMP_PRECISION calc_rms_gsh(COMP_PRECISION *ar, COMP_PRECISION *ai,
			    COMP_PRECISION *br, COMP_PRECISION *bi,int lmax)
{
  int l;
  COMP_PRECISION rms;
  rms=0.0;
  for(l=1;l<=lmax;l++){		/* this is OK, we're summing over zero
				   terms */
    rms += ((COMP_PRECISION)(4*l+2))*
      degree_power_gsh(ar,ai,br,bi,l);
  }
  return(sqrt(rms)/TWO_SQRT_PI);
}
/*
  total power normalized by surface area of sphere, including l =0 term
*/

COMP_PRECISION calc_total_power(COMP_PRECISION *a, COMP_PRECISION *b,int lmax)
{
  int l;
  COMP_PRECISION rms;
  rms=0.0;
  for(l=0;l<=lmax;l++)
    rms += ((COMP_PRECISION)(2*l+1))*degree_power(a,b,l);
  return(sqrt(rms)/TWO_SQRT_PI);
}

/*
  total power normalized by surface area of sphere
*/

COMP_PRECISION calc_total_power_gsh(COMP_PRECISION *ar, COMP_PRECISION *ai,
				    COMP_PRECISION *br, COMP_PRECISION *bi,int lmax)
{
  int l;
  COMP_PRECISION rms;
  rms=0.0;
  for(l=0;l<=lmax;l++)
    rms += ((COMP_PRECISION)(4*l+2))*degree_power_gsh(ar,ai,br,bi,l);
  return(sqrt(rms)/TWO_SQRT_PI);
}

