/* this version of cplknl_anis has merged all existing versions of
 cplknl_anis and can handle the following anisotropic parametrization
 options:

iopt = 0  ACFLN rho
iopt = 1  ph pv eta sv sh rho 
iopt = 2  SV, SH, Piso phi eta rho
iopt = 3  Siso, xi, Piso phi eta rho
iopt = 4  Siso, xi, Piso phi eta rho, 2-psi anisotropy 
iopt = 5  Siso, xi, Piso phi eta rho, 2-psi anisotropy, 4-psi anisotropy

all resulting kernels are for RELATIVE perturbations of the model
parameters BR 08/31/04 

Some bugs have been fixed FM 09/17/04 

Some comments:
- The anisotropic kernels for the 2-psi term (B,H,G) are not computed if s!=0 (ellipticity case) */

#include <math.h>
/*#include "/home/u1/li/prem/prem.h"*/
#include "dimensiony.h"
#include "/data/25/gung/,eigfun/premy_st.h"
/*#include "/home/u1/li/eiglib/struct.h"*/
#include "/data/25/gung/,eigfun/,kernel/,src/eigeny_st.h"
#include "/home/u1/li/iolib/message.h"
/* #include "/scr/04/barbara/anis/struct_anis.h" */
#include "struct_anisall.h"

void cplknl_anisall(s,eig1,eig2,knla,prem,iopt)
int s,iopt;
knla_st_dyn *knla;
prem_st *prem;
eigen_st *eig1,*eig2;
{
static int first=1;
static float *dmy,*res;
float *farray1();
float wcom,w0n2;
int nod,nod1,nod2,npu;
int flag,i;
float L1,L2,S,b1,b2,b11,b12,B1,B2;
float *wk,*g1,*g2,*bb1,*bb2;
int l1,l2;
char q1,q2;
  
  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(cplknl: to be fixed);
  if(q1!=q2) STOP(cplknl: 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;
    nod2=prem->majordisc.nsl;
  }
  else
  { nod1=0;
    nod2=prem->majordisc.nsurf;
  }
  
  L1=(float)((l1+1)*l1);
  L2=(float)((l2+1)*l2);
  S=(float)((s+1)*s);
  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;

  if (s==0) {
     B1=0.5*sqrt(L1*(L1-2));
     B2=0.5*sqrt(L2*(L2-2));
  }

/* define and compute g1 and g2 */

  g1=farray1(nod1,nod2);
  g2=farray1(nod1,nod2);

  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,etan;
    float rka0,rkc0,rkf0,rkn0,rkga0,rkgb0,rkl1,rkn2,rkr1,rkga11,rkga12,
          rkgb11,rkgb12,rkr0;
    rr=prem->radius[nod];
    rr2=rr*rr;
    rho=prem->rho[nod];
    gr=prem->g[nod];
    
    u1=eig1->u[nod];
    v1=eig1->v[nod];
    p1=eig1->ph[nod];
    up1=eig1->up[nod]*rr;
    vp1=eig1->vp[nod]*rr;
    pp1=eig1->php[nod]*rr;
    
    u2=eig2->u[nod];
    v2=eig2->v[nod];
    p2=eig2->ph[nod];
    up2=eig2->up[nod]*rr;
    vp2=eig2->vp[nod]*rr;
    pp2=eig2->php[nod]*rr;
    
    corfac(nod,wcom,q1,prem,&xac,&xf,&xln);
    
    aa=xac*prem->a[nod];
    cc=xac*prem->c[nod];
    ff=xf*prem->f[nod];
    ll=xln*prem->l[nod];
    nn=xln*prem->n[nod];
    
    x1=vp1-v1+u1;
    x2=vp2-v2+u2;
    usq=u1*u2;
    vsq=v1*v2;
    xsq=x1*x2;
   
    f1=2.*u1-L1*v1; if(q1=='T') f1=0.0;
    f2=2.*u2-L2*v2; if(q2=='T') f2=0.0;

    rkga0=rho*usq*S;
    rkgb0=-rho*(f1*u2+f2*u1)*rr;
    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*rr;
    rkgb12=.5*rho*u2*v1*rr;

    g1[nod]=rkga0+b11*rkga11+b12*rkga12;
    g2[nod]=rkgb0+b11*rkgb11+b12*rkgb12;
  }


/* compute the integrals entering in the definition of R(2) 
equation A48 in Woodhouse (1980) BAR05/14/2002 */

  wk=farray1(nod1,nod2);
  dmy=farray1(nod1,nod2);
  res=farray1(nod1,nod2);
  bb1=farray1(nod1,nod2);
  bb2=farray1(nod1,nod2);
  
  flag=0;

  for (nod=nod1;nod<=nod2;nod++)
  { float rr,rr2,xx1,xx2,rx,rx2;

    rr=prem->radius[nod];
    rr2=rr*rr;
    if (q1=='S') {

      for (i=nod;i<=nod2;i++) {
	rx=prem->radius[i];
    	rx2=rx*rx;
/* The normalization factor is not needed, check notes, Federica 5/23/05 */
	wk[i]=((s+1)*g2[i]-rx*g1[i]);
/*	wk[i]=((s+1)*g2[i]-rx*g1[i])*rr2/rx2;*/
/*       if(s) wk[i]*=(pow((double)rr,(double)s))/(pow((double)rx,(double)s));*/
        if(s) {
	  if (rx)/* Changed to avoid division by 0, Federica 5/23/05  */ wk[i]*=(pow((double)rr,(double)s))/(pow((double)rx,(double)s));
	  else wk[i]=0;
	}
      }
      npu=nod2-nod+1;
      pingy_(&flag,&npu,prem->radius+nod,wk+nod,dmy+nod,res+nod);
      bb1[nod]=res[nod2];
      for (i=nod1;i<=nod;i++) {
	rx=prem->radius[i];
/* The normalization factor is not needed, check notes */
	wk[i]=(s*g2[i]+rx*g1[i]);
/*	wk[i]=(s*g2[i]+rx*g1[i])*rr2/rx2;
        wk[i]*=(pow((double)rx,(double)(s+1)))/(pow((double)rr,(double)(s+1))); */
	if (rr) /* Changed to avoid division by 0, Federica 5/23/05  */ wk[i]*=(pow((double)rx,(double)(s+1)))/(pow((double)rr,(double)(s+1)));
	else wk[i]=0;
      }
      npu=nod-nod1+1;
      pingy_(&flag,&npu,prem->radius+nod1,wk+nod1,dmy+nod1,res+nod1);
      bb2[nod]=res[nod];
    } else {
      bb1[nod]=0;bb2[nod]=0.;
    }
  } /* end computation of integrals bb1 and bb2 */

  free_farray1(g1,nod1);
  free_farray1(g2,nod1);
  free_farray1(res,nod1);
  free_farray1(wk,nod1);
  free_farray1(dmy,nod1);

  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,etan;
    float rka0,rkc0,rkf0,rkn0,rkga0,rkgb0,rkl1,rkn2,rkr1,rkga11,rkga12,
          rkgb11,rkgb12,rkr0,rkr22;
    float pv,ph,sv,sh,ksv,ksh;
    float piso,katmp,kctmp;
    float rkh0,rkb0;

    rr=prem->radius[nod];
    rr2=rr*rr;
    rho=prem->rho[nod];
    gr=prem->g[nod];
    
    u1=eig1->u[nod];
    v1=eig1->v[nod];
    p1=eig1->ph[nod];
    up1=eig1->up[nod]*rr;
    vp1=eig1->vp[nod]*rr;
    pp1=eig1->php[nod]*rr;
    
    u2=eig2->u[nod];
    v2=eig2->v[nod];
    p2=eig2->ph[nod];
    up2=eig2->up[nod]*rr;
    vp2=eig2->vp[nod]*rr;
    pp2=eig2->php[nod]*rr;
    
    corfac(nod,wcom,q1,prem,&xac,&xf,&xln);
    
    aa=xac*prem->a[nod];
    cc=xac*prem->c[nod];
    ff=xf*prem->f[nod];
    ll=xln*prem->l[nod];
    nn=xln*prem->n[nod];
    
    x1=vp1-v1+u1;
    x2=vp2-v2+u2;
    usq=u1*u2;
    vsq=v1*v2;
    xsq=x1*x2;
   
    f1=2.*u1-L1*v1; if(q1=='T') f1=0.0;
    f2=2.*u2-L2*v2; if(q2=='T') f2=0.0;
    
    rka0=f1*f2;
    rkc0=up1*up2;
    rkf0=up1*f2+up2*f1;
    rkn0=-rka0;
    rkr0=((8.*rho-w0n2)*usq*rr+pp1*u2+pp2*u1-.5*gr*(4.*usq+f1*u2+f2*u1))*rr;
    rkga0=rho*usq*S;
    rkgb0=-rho*(f1*u2+f2*u1)*rr;
    rkl1=xsq;
    rkn2=vsq;
    rkr1=-w0n2*vsq*rr2+(p1*v2+p2*v1+.5*gr*(u1*v2+u2*v1))*rr;
    rkr22=4./(2.*s+1)*(bb1[nod]-bb2[nod]);
    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*rr;
    rkgb12=.5*rho*u2*v1*rr;
    
    if (iopt>=4 && s==0) {
      rkh0=B2*up1*v2+B1*up2*v1;
      rkb0=B2*f1*v2+B1*f2*v1;
    }
    
    if(q1!=q2)
    { if(q2=='T')
      { rkr1=-rkr1;
        rkl1=-rkl1;
        rkn2=-rkn2;
      }
      else
      { rkga12=-rkga12;
        rkgb12=-rkgb12;
      }
    }
    
    knla->g1[nod]=rkga0+b11*rkga11+b12*rkga12;
    knla->g2[nod]=rkgb0+b11*rkgb11+b12*rkgb12;

    if(q1==q2) knla->Coriolis[nod]=rr2*(vsq+u1*v2+u2*v1);
    else STOP(cplknl: to be fixed 3);

  /* parametrization in ACFLN rho */
  if (iopt==0) { 

    if (q1=='S') {
      float xmu;
      etan=ff/(aa-2.*ll);
      knla->a[nod]=rka0*aa;
      knla->c[nod]=rkc0*cc;
      knla->f[nod]=rkf0*ff;
      xmu=rkn0+b2*rkn2+b1*rkl1+4.*(rka0+rkc0)/3.-(2./3.)*rkf0;
/*    printf(" r %g aa %g rka0 %g cc %g rkc0 %g ff %g rkf0 %g\n",
      rr,aa,rka0,cc,rkc0,ff,rkf0); */
/*    printf(" r %g rkn0 %g b2 %g rkn2 %g b1 %g rkl1 %g xmu %g\n",
      rr,rkn0,b2,rkn2,b1,rkl1,xmu); */
     
    } else {
      etan=0.0;
      knla->a[nod]=0.;
      knla->c[nod]=0.;
      knla->f[nod]=0.;
    }
    
/*  if(nn*ll==0.0) knl->vs[nod]=knl->mu[nod]=0.0; 
    else */

    knla->n[nod]=(rkn0+b2*rkn2)*nn;
    knla->l[nod]=rkl1*b1*ll;
    knla->rho[nod]=(rkr0+b1*rkr1+rkr22)*rho;

  } /* end parametrization in ACFLN rho */

  else if (iopt>=1) {

/****************************    
	if(q1=='S') 
    { float xmu;
	etan=ff/(aa-2.*ll);
	knla->a[nod]=rka0;
	knla->c[nod]=rkc0;
	knla->f[nod]=rkf0;
	xmu=rkn0+b2*rkn2+b1*rkl1+4.*(rka0+rkc0)/3.-(2./3.)*rkf0;
    }
    else 
    {   etan=0.0;
	knla->a[nod]=0.;
	knla->c[nod]=0.;
	knla->f[nod]=0.;
    }			*******************/    

    knla->n[nod]=(rkn0+b2*rkn2);
    knla->l[nod]=rkl1*b1;
    knla->rho[nod]=(rkr0+b1*rkr1+rkr22);

    /*A->ph, C->pv, L->sv, N->sh, rho->rho' p336   gung 07/24/02 */
    /* knla->a multiplied by ph, since it is the kernel
       of dph/ph, same as others */  

    pv=sqrt(cc/rho);
    ph=sqrt(aa/rho);
    sh=sqrt(nn/rho);
    sv=sqrt(ll/rho);
    piso=sqrt((3*cc+8*aa+8*ll+4*ff)/(15.*rho));

    if (q1=='S') {
	etan=ff/(aa-2.*ll);

        knla->rho[nod]=rho*((pv*pv*rkc0)         /* p.336, 9.43 */
	+ph*ph*(rka0+etan*rkf0)
	+(sv*sv*(knla->l[nod]-2*etan*rkf0))
	+sh*sh*(knla->n[nod])
	+knla->rho[nod]);

      katmp=2*rho*ph*(rka0+etan*rkf0);
      kctmp=2*rho*pv*rkc0;
      
      knla->a[nod]=katmp*ph; /*kph*/
      knla->c[nod]=kctmp*pv; /*kpv*/
      knla->l[nod]=sv*2*rho*sv*(knla->l[nod]-2*etan*rkf0); /*ksv*/
      knla->n[nod]=sh*2*rho*sh*(knla->n[nod]); /*ksh*/
      knla->f[nod]=etan*rho*(ph*ph-2*sv*sv)*rkf0; /*keta*/

      if (iopt>=2) {

     /* from K (ph pv eta sv sh rho ) - > K (Piso phi eta sv sh rho) */
     /* kernel a c are replaced by piso and phi 08/2002/ gung */

        knla->a[nod]=pv*kctmp+ph*katmp;
        knla->c[nod]=ph*pv*(4*ph*kctmp-pv*katmp)/(piso*piso*10);
      }

    } else {
     
      etan=0.0;
      knla->rho[nod]= rho*(sv*sv*(knla->l[nod])+sh*sh*(knla->n[nod])	 
	 +knla->rho[nod]);
      knla->l[nod]=sv*2*rho*sv*(knla->l[nod]);
      knla->n[nod]=sh*2*rho*sh*(knla->n[nod]);
	
      knla->a[nod]=0.;
      knla->c[nod]=0.;
      knla->f[nod]=0.;
    }
      
    if (iopt>=3) {

    /* from K (sv sh) - > K (Siso xi) */
    /* kernel l n are replaced by Siso and xi 10/2002/ gung */

      ksv=knla->l[nod];
      ksh=knla->n[nod];
      knla->l[nod]=ksv+ksh;
      knla->n[nod]=(ll*ksh-0.5*nn*ksv)/(2.*ll+nn);
    }
 
  } /* end iopt==1 */

  if (iopt>=4 && s==0) { 
  /* additional 2-psi azimuthal terms B,G,H  
     currently relative values of perturbations:
     - for G with respect to L
     - for B with respect to A
     - for H with respect to F */

    if (q1=='S') {
       knla->g[nod]=b1*rkl1*ll;
       knla->b[nod]=-rkb0*aa;
       knla->h[nod]=-2*rkh0*ff;
    } else {
       knla->g[nod]=-b1*rkl1*ll;
       knla->b[nod]=0;
       knla->h[nod]=0;
    }   
  }

  if (iopt==5) {
  /* additional 4-psi azimuthal term CS */
  /* currently relative values of perturbations with respect to N */
	
    if (q1=='S') knla->cs[nod]=b2*vsq*nn;
    else knla->cs[nod]=-b2*vsq*nn;
  }

 } /* end  loop on nod */  

 free_farray1(bb1,nod1);
 free_farray1(bb2,nod1);
}