/* 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 */

#include <stdio.h>
#include <math.h>
#include "/home/u1/li/iolib/message.h"
#include "dimensiony.h"
#include "/data/25/gung/,eigfun/premy_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 MODETyPE T
#define S 2  
#define OM2 5.2885e-09 /* earth's rotation rate */

float ellcoup(eig1,eig2,prem)
prem_st *prem;
eigen_st *eig1,*eig2;
{
double add,comp,rot;
/*knl_st knl;*/
static int n1,n2,first=1;
static float *integ,*dmy,*res;
float cplknl_disc();
int l1,l2,i,flag,nn,ii;
float dr,el,fac,add1,add2,add3;
float wcom,w0n2;
int nod,nod1,nod2,npu;
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(ell_grav: to be fixed);
  if(q1!=q2) STOP(ell_grav: 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) return((float)(0.0));

  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));

  flag=0;
  add=0.0;
  alphaprime=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;
/*  	printf("nod %d dr %g SMALL %g\n",nod,dr,SMALL); */


    /** 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;
	}
      }
    }
  }

  /* 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);
/*	printf("i %d ell %g rad %g add %g\n",
	i,prem->ell[i],prem->radius[i],add); */

    }
  }

  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);
	/* printf("add %g\n",add); */

  add*=0.666666666666;

  if(eig1->type==3||eig1->type==1){
	alphaprime*=OM2/3.;
   add+=alphaprime; 
	}
  /*add*=fac;*/
  
/*  printf("ell term= %g\n",add); */

  return((float)add);  
}