/* convert from disc.f */

#include <math.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 "/home/u1/li/iolib/message.h"

float cplknl_disc(s,iq,eig1,eig2,prem)
int s,iq;
prem_st *prem;
eigen_st *eig1,*eig2;
{
char q1,q2;
float wcom,w0n2,cst;
int l1,l2,i,idis,iqt;
float L1,L2,S,b1,b2,b11,b12;
float acum[5];
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,ttt;
float rka0,rkc0,rkn0,rkga0,rkgb0,rkl1,rkn2,rkr1,rkr0,rkc11,rkc12,
      rkf11,rkf12;


  wcom=sqrt((double)(eig1->w*eig2->w));
  w0n2=(wcom/prem->normfac.wn); w0n2*=w0n2;
  
  l1=eig1->l;
  l2=eig2->l;
  q1='S'; if(eig1->type==2)  q1='T';
  q2='S'; if(eig2->type==2)  q2='T';
  
  L1=(l1+1)*l1;
  L2=(l2+1)*l2;
  S=(s+1)*s;
  b1=0.5*(L1+L2-S);
  b11=.5*(S+L2-L1);
  b12=.5*(S+L1-L2);
  b2=2.*(b1-1.)*b1-L1*L2;
  
  if(q1!=q2) STOP(cplknl_disc: to be fixed);
  
  for(i=0;i<5;i++) acum[i]=0.0;
  
  rr=prem->radius[iq];
  rr2=rr*rr;
  
  for(idis=1;idis<=3;idis+=2)
  { if(iq!=prem->majordisc.nsurf || idis!=3)
    { iqt=iq+(idis-1)/2;
      gr=prem->g[iqt];
      rho=prem->rho[iqt];
      
      u1=eig1->u[iqt];
      v1=eig1->v[iqt];
      p1=eig1->ph[iqt];
      up1=eig1->up[iqt]*rr;
      vp1=eig1->vp[iqt]*rr;
      pp1=eig1->php[iqt]*rr;
            
      u2=eig2->u[iqt];
      v2=eig2->v[iqt];
      p2=eig2->ph[iqt];
      up2=eig2->up[iqt]*rr;
      vp2=eig2->vp[iqt]*rr;
      pp2=eig2->php[iqt]*rr;

      corfac(iqt,wcom,q1,prem,&xac,&xf,&xln);
    
      aa=xac*prem->a[iqt];
      cc=xac*prem->c[iqt];
      ff=xf*prem->f[iqt];
      ll=xln*prem->l[iqt];
      nn=xln*prem->n[iqt];
    
      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;
      rkn0=-rka0;
      rkr0=((8.*rho-w0n2)*usq*rr+pp1*u2+pp2*u1-.5*gr*(4.*usq+f1*u2+f2*
            u1))*rr;
      rkl1=xsq-vp1*x2-vp2*x1;
      rkn2=vsq;
      rkr1=-w0n2*vsq*rr2+(p1*v2+p2*v1+.5*gr*(u1*v2+u2*v1))*rr;
      
      rkc11=up1*v2;
      rkc12=up2*v1;
      rkf11=f1*v2;
      rkf12=f2*v1;
      
      if(q1!=q2)
      { if(q2=='T')
        { rkn2=-rkn2;
          rkl1=-rkl1;
          rkr1=-rkr1;
        }
        else
        { rkc12=-rkc12;
          rkf12=-rkf12;
        }
      }
      
      ttt=(float)(idis-2);
      
      
      if(q1!=q2) STOP(cplknl_disc: to be fixed);
      
      if(q1=='S')
      { acum[0]+=-ttt*(aa*rka0+cc*rkc0+nn*rkn0+rho*rkr0);
        acum[1]+=-ttt*(ll*rkl1+rho*rkr1);
        acum[2]+=-ttt*(cc*rkc11+ff*rkf11);
        acum[3]+=-ttt*(cc*rkc12+ff*rkf12);
        acum[4]+=-ttt*nn*rkn2;
      }
      else
      { acum[0]+=-ttt*(nn*rkn0+rho*rkr0);
        acum[1]+=-ttt*(ll*rkl1+rho*rkr1);
        acum[2]+=0.0;
        acum[3]+=0.0;
        acum[4]+=-ttt*nn*rkn2;
      }
     
    }
  }
  
/*  Corrected by federica (Jan. 10) according to Mark's email (Jan 5) */
/*  cst=acum[0]+acum[1]*b1+acum[2]*b11+acum[3]*b12+acum[4]*b2;*/
  cst=acum[0]+acum[1]*b1+acum[4]*b2;
  
  return(cst);

}