#include <string.h>
#include <math.h>
#include "../iolib/message.h"
#include "cmplx.h"

float **farray2();
complex *cmplxarray1();

void rotharm(a,lmax,jump,conven,alpha,beta,gamma)
float *a;	/* real harmonics to be transformed */
int lmax;       /* max degree of spherical harmonics */
int jump;       /* block size of each harmonic order */
char *conven;   /* "X^",etc. spherical harmonic convention */
float alpha,beta,gamma;  /* Euler angles */
{
float fac,sign;
float **d;
complex *c2,*c1,sum;
double arg;
int i,t,m,l;
float save;

	if(!strcmp(conven,"X^"))
	{ fac=1.4142136;
	}
	else STOP(rotharm: to be fixed);
		
	for(l=0;l<=lmax;l++)
	{ if(l) 
	  { free_farray2(d,-l+1,l-1,-l+1);
	    free_cmplxarray1(c2,-l+1);
	    free_cmplxarray1(c1,-l+1);
	  }
	  d=farray2(-l,l,-l,l);
	  c2=cmplxarray1(-l,l);
	  c1=cmplxarray1(-l,l);
	  
	  rotmatrix(l,l,beta,d);

	  for(i=0;i<jump;i++)
	  { c1[0].real=a[i+jump*l*l]*fac;
	    c1[0].imag=0.0;
	    sign=1.0;
	    for(t=1;t<=l;t++)
	    { sign=-sign;
	      c1[t].real=a[i+jump*(l*l+2*t-1)];
	      c1[t].imag=a[i+jump*(l*l+2*t)];
	      c1[-t].real=sign*c1[t].real;
	      c1[-t].imag=-sign*c1[t].imag;
	    }
	    
	    for(m=0;m<=l;m++)
	    { sum.real=sum.imag=0.0;
	      for(t=-l;t<=l;t++)
	      { arg=gamma*(float)t;
	        sum.real+=
	          d[t][m]*(cos(arg)*c1[t].real-sin(arg)*c1[t].imag);
	        sum.imag+=
	          d[t][m]*(cos(arg)*c1[t].imag+sin(arg)*c1[t].real);
	      }
	      arg=alpha*(float)m;
	      c2[m].real=cos(arg)*sum.real-sin(arg)*sum.imag;
	      c2[m].imag=cos(arg)*sum.imag+sin(arg)*sum.real;	    
	    }
	    
	    a[i+jump*l*l]=c2[0].real/fac;
	    for(m=1;m<=l;m++)
	    { a[i+jump*(l*l+2*m-1)]=c2[m].real;
	      a[i+jump*(l*l+2*m)]=c2[m].imag;
	    }
	  }
	  
	}
	free_farray2(d,-lmax,lmax,-lmax);
        free_cmplxarray1(c2,-lmax);
        free_cmplxarray1(c1,-lmax);

}