/* Ellipticity coupling parameters,
/* B. Romanowicz 06/21/96  */
/* integration over gravitational potential perturbations and G kernels 
   C.Kuo 4/7/97
	final debugging done 10/13/98 B. Romanowicz (renamed to ellcoup) */ 
/* bug on initialization of up2=0. found by Federica 02/06/05 */
/* bug on initialization of alphaprime=0.0 corrected on 02/07/05 */
/* The computation of the rotational splitting parameter has been added 
   (Federica 02/16/2005)*/

/* The routine ellcoup_rot is the function actually computing the ellipticity
   and rotational splitting parameters. The routine ell is a help function 
   used when calling the C routine ellcoup_rot in a f77 code 
   (e.g. yannos_ST_format800km_3_ellrot.f). This routine is also normalizing
   the ellipticity term in the same way as done for the original XD
   eigenfunctions */
   
#include <stdio.h>
#include <math.h>
#include "/home/u1/li/iolib/message.h"
#include "dimensiony.h"
#include "/data/07/gung/,eigfun/premy_st.h"
/*#include "/data/25/gung/,eigfun/premy_st.h"*/
#include "/data/07/gung/,eigfun/eigeny_st.h"
/*#include "/data/25/gung/,eigfun/,kernel/,src/eigeny_st.h"*/
#include "/data/25/gung/,eigfun/,kernel/,src/kernely_st.h"

#define SMALL 1.e-04
#define S 2  
#define OM2 5.2885e-09 /* earth's rotation rate */

extern float *farray1();
extern void free_farray1();

void ell();
void ellcoup_rot();
void eig_farray1();

void ell_(knot,n,type,l,w,q,cgp,a_vert,a_hori,phis,eigbuf,ellterm,rotterm)
int *knot,*n,*type,*l;
float *w,*q,*cgp,*a_vert,*a_hori,*phis,*ellterm,*rotterm,*eigbuf;
{
eigen_st eig;
prem_st prem;
int i,nlmax,mknots;
float fac,dw2;

  nlmax=NNDLPNT;
  loadpremy_(&nlmax,&(prem.ifanis),&(prem.normfac.rn));
  mknots=prem.majordisc.nsurf;
  if (*knot!=mknots+1) stop("ell: model incompatibility");
  eig_farray1(&eig,mknots);

  eig.n=*n;
  eig.type=*type;
  eig.l=*l;
  eig.w=*w;
  eig.q=*q;
  eig.cgp=*cgp;
  eig.a_vert=*a_vert;
  eig.a_hori=*a_hori;
  eig.phis=*phis;
  
  for (i=0;i<(mknots+1);i++) {
    if (*type==2) {  /* Toroidal modes */
      eig.u[i]=0.;
      eig.up[i]=0.;
      eig.v[i]=eigbuf[i];
      eig.vp[i]=eigbuf[*knot+i];
      eig.ph[i]=0.;
      eig.php[i]=0.;

    } else if (*type==1 || *type==3) {  /* Spheroidal modes */
      eig.u[i]=eigbuf[i];
      eig.up[i]=eigbuf[*knot+i];
      eig.v[i]=eigbuf[2*(*knot)+i];
      eig.vp[i]=eigbuf[3*(*knot)+i];
      eig.ph[i]=eigbuf[4*(*knot)+i];
      eig.php[i]=eigbuf[5*(*knot)+i];
    }
  }

  if (*l) {
    ellcoup_rot(&eig,&eig,&prem,ellterm,rotterm);
    /* Normalization as in the original XD eigenfunctions */
    dw2=eig.w*eig.w*(*ellterm);
    fac=(float)((2*eig.l+3)*(2*eig.l-1))/(float)(eig.l*(eig.l+1));
    fac*=2.*eig.w*eig.w;
    fac=1./fac;
    *ellterm=fac*dw2;
  } else {
    *ellterm=0.;
    *rotterm=0.;
  }
  
  free_farray1(eig.u,0);
  free_farray1(eig.up,0);
  free_farray1(eig.v,0);
  free_farray1(eig.vp,0);
  free_farray1(eig.ph,0);
  free_farray1(eig.php,0);
  
  return;
}   
       
/**************************************************************************/
/* Computes and returns the ellipticity (ellterm) and rotational (rotterm)
   splitting parameters. To obtain the same ellipticity term as the one stored
   in the original XD eigenfuctions (e.g. prem.ell) and stored in the structure
   sprem_st (/data/07/gung/,eigfun/premy_st.h), the value of ellterm return by
   this routine needs to be multiplied by some normalization factors (see
   below), as done e.g. in the help routine ell.
   The value for the rotational slitting parameter returned by this
   subroutine (rotterm) has already the same normalization as the one stored
   in the original XD eigenfunctions.

   Used normalization:
   -------------------
   float ellterm,rotterm;
   sprem_st prem;
   prem_st prem_ell;
   eigen_st eigen;

   if (l) {
     ellcoup_rot(&eigen,&eigen,&prem_ell,&ellterm,&rotterm);
     dw2=eigen.w*eigen.w*ellterm; 
     fac=(float)((2*l+3)*(2*l-1))/(float)(l*(l+1));
     fac*=2.*eigen.w*eigen.w;
     fac=1./fac;
     prem.ell=fac*dw2;
     prem.rot=rotterm;
   } else {
     prem.ell=0.;
     prem.rot=0.;
   }
*/
      
void ellcoup_rot(eig1,eig2,prem,ellterm,rotterm)
prem_st *prem;
eigen_st *eig1,*eig2;
float *ellterm,*rotterm;
{
double add,beta,beta1,beta2,rot;
float cplknl_disc();
int l1,l2,i,ii;
float dr,el,fac,add2;
float wcom,w0n2;
int nod,nod1,nod2;
float L1,L2,s,b1,b2,b11,b12;
char q1,q2;

/**** variables for precision ****/
static double wt[5]={0.236926885056189, 0.478628670499366,
              0.568888888888889, 0.478628670499366,0.236926885056189};
static double xi[5]={-0.90617984598664,-0.538469310105683,
              0.000000000000000, 0.538469310105683,0.906179845938664};
float r,halfdr,hr,hsq,hcu,delldr,delldr1;
float Aell,Bell,Av1,Bv1,Av2,Bv2,Au1,Bu1,Au2,Bu2,Ap1,Bp1,Ap2,Bp2;
int subnod;
/*********************************/
float eta,phi,dphidr,alphaprime;

  l1=eig1->l;
  l2=eig2->l;
  q1='S'; if (eig1->type==2) q1='T';
  q2='S'; if (eig1->type==2) q2='T';

  if (!eig1->type || !eig2->type) STOP(ellcoup_rot: to be fixed);
  if (q1!=q2) STOP(ellcoup_rot: to be fixed 2);
  
  wcom=sqrt((double)(eig1->w*eig2->w));
  w0n2=(wcom/prem->normfac.wn);
  w0n2*=w0n2;
  
  if (q1=='T') {
/*    nod1=prem->majordisc.noc; */
    nod1=prem->majordisc.noc+1;
    nod2=prem->majordisc.nsl; 
  } else {
    nod1=1;
    nod2=prem->majordisc.nsurf;
  }

  /* Ellipticity selection rules */
  if (l1!=l2 && l1-2!=l2 && l1+2!=l2) {
    *ellterm=0.0;
    *rotterm=0.0;
    goto finished1;
  }

  L1=(float)((l1+1)*l1);
  L2=(float)((l2+1)*l2);
  s=6.;		/* s=S(S+1).  S=2  */
  b1=0.5*(L1+L2-s);
  b11=0.5*(s+L2-L1);
  b12=0.5*(s+L1-L2);
  b2=2.*(b1-1.)*b1-L1*L2;
  fac=0.5*L1/(float)((2*l1+3))/(float)((2*l1-1));

  add=0.0;
  alphaprime=0.0;
  beta=0.0;
  beta1=0.0;
  beta2=0.0;

  for (nod=nod1;nod<nod2;nod++) { 
  float rr,rr2,rho,gr,u1,v1,up1,vp1,p1,pp1,u2,v2,up2,vp2,p2,pp2;
  float xac,xf,xln,aa,cc,nn,ll,ff,x1,x2,usq,vsq,xsq,f1,f2;
  float rka0,rkc0,rkf0,rkn0,rkl1,rkn2,rkr1,rkr0,rkr2,rkr11,rkr22;
  float rkga0,rkgb0,rkga11,rkga12,rkgb11,rkgb1,rkgb12;
  float dLdr,dAdr,dNdr,dFdr,drhodr,dCdr;
  float xx1,xx2,rx,rx2;
  float dgrdr,dgrdr1,Agr,Bgr,G1,G2;
    rr=prem->radius[nod];
    rr2=rr*rr;
    dr=prem->radius[nod+1]-rr;

    /** start XD interpolation section **/
    if (dr>SMALL) {
      halfdr=0.5*dr;
      hr=1./dr;
      hsq=hr*hr;
      hcu=hr*hsq;

      rho=prem->rho[nod];
      gr=prem->g[nod];
      eta=prem->eta[nod];

      dgrdr= 4.*prem->rho[nod]-2.*prem->g[nod]/prem->radius[nod];
      dgrdr1=4.*prem->rho[nod+1]-2.*prem->g[nod+1]/prem->radius[nod+1];

      if (nod) delldr=prem->ell[nod]*prem->eta[nod]/prem->radius[nod];
      delldr1=prem->ell[nod+1]*prem->eta[nod+1]/prem->radius[nod+1];
   
      /** coefficients for precision **/
      Agr=(dgrdr+dgrdr1)*hsq+2.*(prem->g[nod]-prem->g[nod+1])*hcu;
      Bgr=-(2.*dgrdr+dgrdr1)*hr-3.*(prem->g[nod]-prem->g[nod+1])*hsq;

      Aell=(delldr+delldr1)*hsq+
           2.*(prem->ell[nod]-prem->ell[nod+1])*hcu;
      Bell=-(2.*delldr+delldr1)*hr
           -3.*(prem->ell[nod]-prem->ell[nod+1])*hsq;
      Av1=(eig1->vp[nod]+eig1->vp[nod+1])*hsq+
            2.*(eig1->v[nod]-eig1->v[nod+1])*hcu;
      Bv1=-(2.*eig1->vp[nod]+eig1->vp[nod+1])*hr
          -3.*(eig1->v[nod]-eig1->v[nod+1])*hsq;
      Av2=(eig2->vp[nod]+eig2->vp[nod+1])*hsq+
            2.*(eig2->v[nod]-eig2->v[nod+1])*hcu;
      Bv2=-(2.*eig2->vp[nod]+eig2->vp[nod+1])*hr
          -3.*(eig2->v[nod]-eig2->v[nod+1])*hsq;

      if (eig1->type==3||eig1->type==1) {
      Au1=(eig1->up[nod]+eig1->up[nod+1])*hsq+
            2.*(eig1->u[nod]-eig1->u[nod+1])*hcu;
      Bu1=-(2.*eig1->up[nod]+eig1->up[nod+1])*hr
          -3.*(eig1->u[nod]-eig1->u[nod+1])*hsq;
      Au2=(eig2->up[nod]+eig2->up[nod+1])*hsq+
            2.*(eig2->u[nod]-eig2->u[nod+1])*hcu;
      Bu2=-(2.*eig2->up[nod]+eig2->up[nod+1])*hr
          -3.*(eig2->u[nod]-eig2->u[nod+1])*hsq;
      Ap1=(eig1->php[nod]+eig1->php[nod+1])*hsq+
            2.*(eig1->ph[nod]-eig1->ph[nod+1])*hcu;
      Bp1=-(2.*eig1->php[nod]+eig1->php[nod+1])*hr
          -3.*(eig1->ph[nod]-eig1->ph[nod+1])*hsq;
      Ap2=(eig2->php[nod]+eig2->php[nod+1])*hsq+
            2.*(eig2->ph[nod]-eig2->ph[nod+1])*hcu;
      Bp2=-(2.*eig2->php[nod]+eig2->php[nod+1])*hr
          -3.*(eig2->ph[nod]-eig2->ph[nod+1])*hsq;
      }
 
      corfac(nod,wcom,q1,prem,&xac,&xf,&xln);

      /** interpolate over 5 pts for each node **/
      for (subnod=0;subnod<5;subnod++) { 
      float t,r2;
        t=halfdr*(xi[subnod]+1.);
        r=rr+t;
	r2=r*r;

        el=prem->ell[nod]+t*(delldr+t*(Bell+t*Aell));

        v1=eig1->v[nod]+t*(eig1->vp[nod]+t*(Bv1+t*Av1));
        vp1=(eig1->vp[nod]+t*(2.*Bv1+3.*t*Av1))*r;
        v2=eig2->v[nod]+t*(eig2->vp[nod]+t*(Bv2+t*Av2));
        vp2=(eig2->vp[nod]+t*(2.*Bv2+3.*t*Av2))*r;

        dLdr=xln*(prem->ql[nod][0]+
        t*(2.*prem->ql[nod][1]+3.*t*prem->ql[nod][2]));
        dNdr=xln*(prem->qn[nod][0]+
        t*(2.*prem->qn[nod][1]+3.*t*prem->qn[nod][2]));
        drhodr=prem->qrho[nod][0]+
        t*(2.*prem->qrho[nod][1]+3.*t*prem->qrho[nod][2]);

	dAdr=0.;
	dCdr=0.;
	dFdr=0.;
	u1=0.;up1=0.;u2=0.;up2=0.;p1=0.;pp1=0.;p2=0.;pp2=0.;  

	if (eig1->type==3 || eig1->type==1) {

        u1=eig1->u[nod]+t*(eig1->up[nod]+t*(Bu1+t*Au1));
        up1=(eig1->up[nod]+t*(2.*Bu1+3.*t*Au1))*r;
        p1=eig1->ph[nod]+t*(eig1->php[nod]+t*(Bp1+t*Ap1));
        pp1=(eig1->php[nod]+t*(2.*Bp1+3.*t*Ap1))*r;
        u2=eig2->u[nod]+t*(eig2->up[nod]+t*(Bu2+t*Au2));
        up2=(eig2->up[nod]+t*(2.*Bu2+3.*t*Au2))*r;
        p2=eig2->ph[nod]+t*(eig2->php[nod]+t*(Bp2+t*Ap2));
        pp2=(eig2->php[nod]+t*(2.*Bp2+3.*t*Ap2))*r;
        dAdr=xac*(prem->qa[nod][0]+
        t*(2.*prem->qa[nod][1]+3.*t*prem->qa[nod][2]));
        dCdr=xac*(prem->qc[nod][0]+
        t*(2.*prem->qc[nod][1]+3.*t*prem->qc[nod][2]));
        dFdr=xf*(prem->qf[nod][0]+
        t*(2.*prem->qf[nod][1]+3.*t*prem->qf[nod][2]));
	}
 
	x1=vp1-v1+u1;
       	x2=vp2-v2+u2;
        vsq=v1*v2;
        xsq=x1*x2;
        usq=u1*u2;
    
        f1=2.*u1-L1*v1; if (q1=='T') f1=0.0;
        f2=2.*u2-L2*v2; if (q2=='T') f2=0.0;
 
/*BR: replace rr by r and rr2 by r2 in what follows to match ellcoup */
   
        rka0=f1*f2;
        rkc0=up1*up2;
        rkf0=up1*f2+up2*f1;
        rkn0=-rka0;
        rkr0=((8.*rho-w0n2)*usq*r+pp1*u2+pp2*u1-.5*gr*(4.*usq+f1*u2+f2*u1))*r;
        rkl1=xsq*b1;
        rkn2=vsq*b2;
        rkr1=-w0n2*vsq*r2*b1+(p1*v2+p2*v1+.5*gr*(u1*v2+u2*v1))*b1*r;

        rkr2=rkr0+rkr1;
	
        add+=wt[subnod]*(dAdr*rka0+dCdr*rkc0+dFdr*rkf0+dNdr*rkn0)*r*el*halfdr;
        add+=wt[subnod]*(dLdr*rkl1+dNdr*rkn2)*r*el*halfdr;
        add+=wt[subnod]*(drhodr*rkr2)*r*el*halfdr;

/*	if(nod==nod1+5)printf("rlk1 %g rkn2 %g rkn0 %g add %g %g %g %g %g %g %g\n",
	rkl1,rkn2,rkn0,drhodr,rkr2,r,el,halfdr,wt[subnod],add); */

	if (eig1->type==3 || eig1->type==1) {

        phi=rr*el*gr-0.5*r2*OM2; 
        dphidr=(gr*(eta-1.)+4.*rho*r)*el-OM2*r;
        rkga0=rho*usq*s;
        rkgb0=-rho*(f1*u2+f2*u1)*r;
        rkga11=.5*rho*(u1*vp2+(u1-up1-2.*f1)*v2);
        rkga12=.5*rho*(u2*vp1+(u2-up2-2.*f2)*v1);
        rkgb11=.5*rho*u1*v2*r;
        rkgb12=.5*rho*u2*v1*r;
        G1=rkga0+b11*rkga11+b12*rkga12;
        G2=rkgb0+b11*rkgb11+b12*rkgb12;
 
        add+=wt[subnod]*G1*phi*halfdr;
        add+=wt[subnod]*G2*dphidr*halfdr;
 
        alphaprime+=wt[subnod]*G1*r2*halfdr;
        alphaprime+=wt[subnod]*G2*2.*r*halfdr;
	
        /* Beta for spheroidal modes p.362 Dahlen 1968 */
	/* This is basically like eq. 98 of WD78 but with a 
	different normalization */
	beta1+=wt[subnod]*rho*(vsq+u1*v2+u2*v1)*r2*halfdr;
	beta2+=wt[subnod]*rho*(usq+L1*vsq)*r2*halfdr;
	}
      }
    }
  }

  /* Rotational splitting parameter
     The value stored in prem.rot (here rotterm), compared to WD, is normalized by w */
     
  /* Spheroidal modes (p.362 Dahlen 1968 and WD eq.97) */ 
       
  if (eig1->type==1 || eig1->type==3) {
     beta=beta1/beta2; 
     rot=beta*sqrt(OM2)/eig1->w; 
  }
  
  /* Toroidal modes (WD eq.89) */ 

  if (eig1->type==2) { 
     beta=1./(l1*(l1+1)); /* (see top pag. 350 WD) */
     rot=beta*sqrt(OM2)/eig1->w;
  }

  *rotterm=((float)(rot));
     
  /* add discontinuity contributions */

  for (i=nod1-1;i<nod2-1;i++) {  
    if (prem->radius[i+1]-prem->radius[i]<SMALL) { 
/*  float acum,ttt;
    int idis,iqt;
       acum=0.; */
       add-=prem->ell[i]*prem->radius[i]*cplknl_disc(S,i,eig1,eig2,prem);
    }
  }
  
  if (q1=='T') ii=prem->majordisc.nsl;
  if (q1=='S') ii=prem->majordisc.nsurf;

/*  add2=cplknl_disc(2,ii,eig1,eig2,prem); */

  add-=prem->ell[ii]*prem->radius[ii]*cplknl_disc(2,ii,eig1,eig2,prem);
  add*=0.666666666666;

  if (eig1->type==3||eig1->type==1) {
    alphaprime*=OM2/3.;
    add+=alphaprime; 
  }
  /*add*=fac;*/
    
  *ellterm=((float)(add));
  finished1:; 
}

/************************************************/
void eig_farray1(eign,nts)
     eigen_st *eign;
     int nts;
{
  float *farray1();   
  eign->v=farray1(0,nts);
  eign->vp=farray1(0,nts);
  eign->u=farray1(0,nts);
  eign->up=farray1(0,nts);
  eign->ph=farray1(0,nts);
  eign->php=farray1(0,nts);
}