#include "CalcMethod.h"
#include "LEquations.h"
#include "Matrix1.h"

CalcMethod::CalcMethod(void)
{
}

CalcMethod::~CalcMethod(void)
{
}

//\/\/\/\/\/\/\/\/\/\/\/\/\/\/\\/\/\/\/\/\\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\/\\/\/\/\/\/\/\/\

float	PI2=8.*atan(1.);
float	PI7=1./sqrt(2.);
void CalcMethod::stl_dt(float t,float de,int ifmin,int ifmax,float dfinterval,float rb,float R_TOOL,float DEN_BORE_FLUID, float VIS, float POR,float dkb,float dkbf,float PERM,float vbf,float VS,float FMST,float SPACEE,complex<float> kew[512],float FWFSYN[1024],complex<float>CFDWV0[1024],float dtf[3])
{
    int i;
    float s1,s2,s3,s4,ff,ww,QST,P,kdd,PI=3.1415926;
    complex<float>ai,ar,K,ftemp,PERMw,fe,ke,D,kee,BC,zx,wi[1024],wk[1024];
    s1=0.0;s2=0.0;s3=s4=0.0;
    ai=complex<float>(0.0,1.0);
    ar=complex<float>(1.0,0.0);
    //		 AA=complex<float>((8*3*PERM/POR),0.);

    for(i =ifmin;i<ifmax;i++)
    {
        ff=i*dfinterval;
        ww=2*PI*ff;
        ke=kew[i];
        ftemp=complex<float>(t*PERM*DEN_BORE_FLUID*ww/(VIS*POR),0.0)*ai;
        kee=sqrt(ar-ftemp/complex<float>(2.,0.0))-ftemp;
        PERMw=PERM/kee;//动态渗透率
        D=PERMw*dkbf/(VIS*POR*(de));//流体扩散系数
        if(real(ke)-ww/vbf>0)
        {
            fe=sqrt(ke*ke-ww*ww/vbf/vbf);//径向波数
        }
        else
        {
            fe=sqrt(-ke*ke+ww*ww/vbf/vbf);
        }
        CFB(fe*rb,wi,wk,0);
        BC=fe*rb*wi[1]/wi[0];
        zx=rb*sqrt(complex<float>(-1.*ww,0.0)*ai/D+ke*ke);
        CFB(zx,wi,wk,0);

        ftemp=complex<float>(2.*DEN_BORE_FLUID*ww/((rb*rb-R_TOOL*R_TOOL)*VIS),0.0)*ai*PERMw
            *zx/(complex<float>(1.,0.0)+pow(BC,complex<float>(VS/vbf,0.0)))
            *wk[1]/wk[0];
        //       if(real(ftemp)>0.0)ftemp=complex<float>(-real(ftemp),imag(ftemp));
        K=sqrt(ke*ke+ftemp); //波数

        QST=0.5*real(K)/imag(K);
        K=K*(ar+1/QST*ai/(complex<float>(2.)))*(ar+1/QST*complex<float>(log(ff/FMST))/(PI));

        P=real(ke)/sqrt(real(K)*real(K)+imag(K)*imag(K));

        kdd=real(K)*SPACEE;//qrt(real(ai*K*complex<float>(SPACEE,0.0))*real(ai*K*complex<float>(SPACEE,0.0))+imag(ai*K*complex<float>(SPACEE,0.0))*imag(ai*K*complex<float>(SPACEE,0.0)));//real(exp(ai*K*SPACE));
        //			 aa[i]=fabs(real(K/ke)*FWFSYN[i]*real(CFDWV0[i])*real(exp(ai*K*SPACEE)));

        s1=P*FWFSYN[i]*kdd+s1;
        s2=s2+i*P*FWFSYN[i]*kdd;
        s3=s3+ww*ww*FWFSYN[i]*FWFSYN[i];//real(CFDWV0[i]);
        s4=s4+(kdd/ww-real(ke)*SPACEE/ww)*ww*ww*FWFSYN[i]*FWFSYN[i];//real(CFDWV0[i]);

    }
    dtf[0]=s2/s1*dfinterval;
    dtf[1]=s4/s3*1000000.;
    dtf[2]=QST;

}

float CalcMethod::prr_N(float wave_R[8][1024],float DT_INTERVAL,float TIME_L,float SPACE_RLEV,float time,float dt,int MMM)
{
    double  s1,s2,s3,s4;
    int iflag;
    float t,awave,aij,rr1;
    int ij,jj;
    s2=0.0;
    s4=0.0;
    iflag=0;
    if(fabs(time-dt*12)<1600)
    {
        for(t=0;t<TIME_L;t=t+DT_INTERVAL)
        {
            s1=0.0;s3=0.0;
            for(ij=0;ij<MMM;ij++)
            {
                aij=(t+time+ij*dt*SPACE_RLEV)/DT_INTERVAL;

                jj=aij;
                if(jj+1>660)iflag=1;
                if(aij>float(jj))
                {
                    awave=(1-aij+jj)*wave_R[ij][jj]+(aij-jj)*wave_R[ij][jj+1];
                }
                else
                {
                    awave=(aij-jj+1)*wave_R[ij][jj]+(jj-aij)*wave_R[ij][jj-1];

                }
                s1=s1+sqrt(sqrt(fabs(awave)))*fabs(awave)/awave;
                s3=s3+sqrt(sqrt(fabs(awave)));
            }
            s2=s2+s1*s1*s1*s1;
            s4=s4+s3*s3*s3*s3;
        }

        rr1=sqrt(sqrt(s2/s4));
    }
    else
    {
        rr1=0.0;
    }
    if(iflag==1)rr1=0.0;
    return rr1;
}

//
//
//以下为子程序
///
////
///
complex<float> CalcMethod::DET(float vp, float vs,float vf,float den,float R, float w ,int IONLY)
{
    return 0;
}
void CalcMethod::CFB(complex<float> Z, complex<float> wi[2], complex<float> wk[2],  int IONLY)
{
    /*
    SUBROUTINE CFB(Z,I0,I1,K0,K1,IONLY)

    MW114420
    C      SUBROUTINE TO CALCULATE MODIFIED BESSEL FUNCTIONS                MW114430
    C      I0,I1,K0,AND K1 FOR COMPLEX ARGUMENTS                            MW114440
    C      EQUATION NUMBERS REFER TO EQUATIONS IN                           MW114450
    C      ABRAMOWITZ AND STEGUN.                                           MW114460
    C      IONLY IS SET TO 1 IF ONE WANTS TO CALCULATE I0 AND I1 ONLY.      MW114470
    */
    complex<float>   Z2O4,ZT8,TOP,TERM,SUM;
    complex<float>   BF[100],I0,I1,K0,K1;
    complex<float>    temp1,temp2;
    float GAMMA,PIO2,TWOPI,PIT3O2,ORDER,SIGN,AMU,FACTK,SUMINV,X,Y,B;
    float CONV ;
    //	float yyi,yyr;
    GAMMA=.5772156649015328;
    PIT3O2=4.712388981;
    PIO2=1.570796326794896;
    TWOPI=6.283185307179586;
    CONV=1.0e-12;
    //C      PI=3.141592653589793
    int I,II,j,JJ,jj,N;
    //C      CALCULATE I0 AND I1 USING RECURRENCE RELATIONSHIP EQN 9.6.26...
    if (fabs(real(Z))<=15.0)
    {
        N=int(8.500+1.500*fabs(real(Z))) ;
        BF[N+2]=complex<float> (0.00,0.00);
        BF[N+1]=complex<float> (1.0e-20,1.0e-20);

        temp2=complex<float> (2.0,0.0);

        SUM=complex<float> (0.00,0.00) ;

        for(I=N;I>0;I--)
        {
            X=float(I);
            temp1=complex<float> (X);
            BF[I]=temp1*temp2/Z*BF[I+1]+BF[I+2];
            SUM=SUM+BF[I];
        }

        //C      NORMALIZE USING EQN 9.6.37
        //MW114710
        SUM=exp(Z)/(BF[1]+temp2*(SUM-BF[1])) ;
        I0=SUM*BF[1]  ;
        I1=SUM*BF[2] ;

        //C      ... OR AYMPTOTIC SERIES EQN 9.7.1
    }
    else
    {
        temp1=complex<float> (8.0,0.0);
        ZT8=Z*temp1 ;
        for(I=1;I<=2 ;I++)
        {
            ORDER=float(I-1) ;
            SIGN=-1.00;
            AMU=4.00*ORDER*ORDER ;
            SUM=complex<float> (1.00,0.00)  ;
            TERM=complex<float> (1.00,0.00) ;
            Y=1.00    ;

            for( j=1;j<=15;j=j+2)
            {
                X=float(j) ;
                TERM=TERM* (complex<float>)(AMU-X*X)/ZT8/(complex<float>)(Y );

                if (fabs(real(TERM))<CONV) goto L_30;
                SUM=SUM+TERM*SIGN;

                Y=Y+1.00 ;
                SIGN=SIGN*(-1.00);
            }

L_30:          if (ORDER==0)
               {

                   I0=exp(Z)/sqrt(TWOPI*Z)*SUM ;

               }
               else
               {
                   I1=exp(Z)/sqrt(TWOPI*Z)*SUM ;

               }

        }
    }
    wi[0]=I0;
    wi[1]=I1;
    //C      CHECK TO SEE IF ONLY I'S ARE NEEDED
    if(IONLY==1)
    {

        return ;
    }

    //C      CALCULATE K0 USING SERIES...                                     MW115090
    //C      ASCENDING SERIES EQN 9.6.13...                                   MW115100

    temp2=complex<float> (2.0,0.0);
    if(fabs(real(Z))<5.0)
    {
        K0=(log(Z/temp2)+complex<float>(GAMMA,0.0))*complex<float>(-1.,0.0)*I0;
        Z2O4=Z*Z/complex<float>(4.00,0.00) ;
        FACTK=1.00;
        TOP=complex<float> (1.00,0.00);
        SUMINV=0.0 ;
        for(jj=1;jj<25;jj++)
        {
            for( JJ=1;JJ<50;JJ++)
            {
                B=float(JJ);
                FACTK=FACTK*B;
                if(JJ>20)
                {
                    float nnn=0;
                }
                TOP=TOP*Z2O4;
                SUMINV=SUMINV+1./B ;
                TERM=TOP/ complex<float>(FACTK,0.0)/complex<float>(FACTK,0.0)*complex<float>(SUMINV,0.0);
                K0=K0+TERM ;
                if (fabs(real(TERM))<CONV) goto L_80;
            }
        }

        //      ...OR ASYMPTOTIC SERIES EQN 9.7.2
    }
    else
    {
        ZT8=Z*complex<float>(8.,0.) ;
        SUM=complex<float> (1.0,0.0) ;
        TERM=complex<float> (1.0,0.0);
        Y=1.0  ;

        for( II=1;II<30;II=II+2)
        {
            X=float(II) ;
            TERM=-TERM*complex<float>(X,0.0)*complex<float>(X,0.0)/ZT8/complex<float>(Y,0.0);
            SUM=SUM+TERM;
            if (fabs(real(TERM))<CONV) goto L_70 ;
            Y=Y+1.00 ;
        }

L_70:    K0=sqrt(PIO2/Z)/exp(Z)*SUM;

    }
    //      FIND K1 USING THE WRONSKIAN EQN 9.6.15
L_80:
    K1=(complex<float>(1.00,0.0)/Z-I1*K0)/I0 ;
    wk[0]=K0;
    wk[1]=K1;
}

///
int CalcMethod::amp(float x[1024], int nl, float xr[1024],int i0)
{
    int i,nfft,mm;
    float xi[1024],xmax;
    CString ss;
    nfft=1024;
    for(i=0;i<1024;i++)xi[i]=0.;
    for(i=0;i<nl;i++)
    {
        xr[i]=x[i];
        xi[i]=0.f;
    }
    for(i=nl;i<nfft;i++)
    {
        xr[i]=0.f;
        xi[i]=0.f;
    }

    if(i0==0)
    {
        fft842(i0,nfft,xr,xi);

    }
    if(i0==1)
    {
        fft842(0,nfft,xr,xi);
        for(i=0;i<nfft;i++)
        {
            float a=sqrt(xr[i]*xr[i]+xi[i]*xi[i]);
            float b=atan(xi[i]/xr[i]);
            xr[i]=a;
            xi[i]=b;
        }
        fft842(1,nfft,xr,xi);
    }
    mm=0;
    if(i0==0)
    {
        for(i=0;i<nfft;i++)
            xr[i]=sqrt(xr[i]*xr[i]+xi[i]*xi[i])/nfft;
        xmax=xr[0];

        for(i=0;i<(nfft/2);i++)
        {
            if(xmax < xr[i])
         {
             xmax=xr[i];
             mm=i;
         }
        }
    }
    return mm;
}

void CalcMethod::fft842(int INN, int N, float *X, float *Y)
{
    int I,IJ,JI,IPASS,IW,J,J1,J2,J3,J4,J5,J6,J7,J8,J9,J10,J11,J12,J13,J14,
        LENGT,M,M1,M2,M3,M4,N2POW,N8POW,NT,NTHPO,NXTLT,L[15];
    float FI,FN,R;
    float XX[1025],YY[1025];
    IW=0;
    //	PI2=8.*atan(1.);
    //	PI7=1./sqrt(2.);
    for(I=1;I<=15;I++)
    {
        M=I;
        NT=pow(2.0,I);
        if(N == NT)	goto L_20;
    }
    //	exit(0);
L_20:
    N2POW=M;
    NTHPO=N;
    FN=1.0f*NTHPO;
    if(INN != 1)
    {
        for(I=0;I<NTHPO;I++)
            Y[I]=-Y[I];
    }
    N8POW=N2POW/3;
    if(N8POW != 0)
    {
        for(IPASS=0;IPASS<N8POW;IPASS++)
        {
            NXTLT=pow(2.0,N2POW-3*(IPASS+1));
            LENGT=8*NXTLT;
            r8tx(NXTLT,NTHPO,LENGT,X,Y);//
            //X+NXTLT,X+NXTLT*2,X+NXTLT*3,X+NXTLT*4,X+NXTLT*5,X+NXTLT*6,
            //	X+NXTLT*7,Y,Y+NXTLT,Y+NXTLT*2,Y+NXTLT*3,Y+NXTLT*4,Y+NXTLT*5,Y+NXTLT*6,Y+NXTLT*7);
        }
    }
    M1=N2POW-3*N8POW-1;
    if(M1==0)
        r2tx(NTHPO,X,Y);//X+1,Y,Y+1);
    else if(M1>0)
        r4tx(NTHPO,X,Y);//X+1,X+2,X+3,Y,Y+1,Y+2,Y+3);
    for(J=1;J<=15;J++)
    {
        L[J-1]=1;
        M2=J-N2POW;
        if(M2 <=0 )
        {
            M3=N2POW-J+1;
            L[J-1]=pow(2.0,M3);
        }
    }
    IJ=1;
    for(J1=1;J1<=N;J1++)
    {
        XX[J1]=X[J1-1];
        YY[J1]=Y[J1-1];
    }
    for(J1=1;J1<=L[14];J1++)
    {
        for(J2=J1;J2<=L[13];J2+=L[14])
        {
            for(J3=J2;J3<=L[12];J3+=L[13])
            {
                for(J4=J3;J4<=L[11];J4+=L[12])
                {
                    for(J5=J4;J5<=L[10];J5+=L[11])
                    {
                        for(J6=J5;J6<=L[9];J6+=L[10])
                        {
                            for(J7=J6;J7<=L[8];J7+=L[9])
                            {
                                for(J8=J7;J8<=L[7];J8+=L[8])
                                {
                                    for(J9=J8;J9<=L[6];J9+=L[7])
                                    {
                                        for(J10=J9;J10<=L[5];J10+=L[6])
                                        {
                                            for(J11=J10;J11<=L[4];J11+=L[5])
                                            {
                                                for(J12=J11;J12<=L[3];J12+=L[4])
                                                {
                                                    for(J13=J12;J13<=L[2];J13+=L[3])
                                                    {
                                                        for(J14=J13;J14<=L[1];J14+=L[2])
                                                        {
                                                            for(JI=J14;JI<=L[0];JI+=L[1])
                                                            {
                                                                M4=IJ-JI;
                                                                if(M4<0)
                                                                {
                                                                    R=XX[IJ];
                                                                    XX[IJ]=XX[JI];
                                                                    XX[JI]=R;
                                                                    FI=YY[IJ];
                                                                    YY[IJ]=YY[JI];
                                                                    YY[JI]=FI;
                                                                }
                                                                IJ++;
                                                            }
                                                        }
                                                    }
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
    //
    for(J1=1;J1<=N;J1++)
    {
        X[J1-1]=XX[J1];
        Y[J1-1]=YY[J1];
    }
    if(INN == 1)
    {
        for(I=0;I<NTHPO;I++)
        {
            X[I]=X[I]/FN;
            Y[I]=Y[I]/FN;
        }
    }
    else
    {
        for(I=0;I<NTHPO;I++)
            Y[I]=-Y[I];
    }
}

void CalcMethod::r2tx(int NTHPO, float *CR0, float *CI0)//1, float *CI0, float *CI1)
{
    int K;
    float FI1,R1;

    for(K=0;K<NTHPO;K+=2)
    {
        R1=CR0[K]+CR0[K+1];
        CR0[K+1]=CR0[K]-CR0[K+1];
        CR0[K]=R1;
        FI1=CI0[K]+CI0[K+1];
        CI0[K+1]=CI0[K]-CI0[K+1];
        CI0[K]=FI1;
    }
    return;
}

void CalcMethod::r4tx(int NTHPO, float *CR0, float *CI0)//, float *CR2, float *CR3, float *CI0, float *CI1, float *CI2, float *CI3)
{
    int K;
    float FI1,FI2,FI3,FI4,R1,R2,R3,R4;

    for(K=0;K<NTHPO;K+=4)
    {
        R1=CR0[K]+CR0[K+2];
        R2=CR0[K]-CR0[K+2];
        R3=CR0[K+1]+CR0[K+3];
        R4=CR0[K+1]-CR0[K+3];
        FI1=CI0[K]+CI0[K+2];
        FI2=CI0[K]-CI0[K+2];
        FI3=CI0[K+1]+CI0[K+3];
        FI4=CI0[K+1]-CI0[K+3];
        CR0[K]=R1+R3;
        CI0[K]=FI1+FI3;
        CR0[K+1]=R1-R3;
        CI0[K+1]=FI1-FI3;
        CR0[K+2]=R2-FI4;
        CI0[K+2]=FI2+R4;
        CR0[K+3]=R2+FI4;
        CI0[K+3]=FI2-R4;
    }
}
void CalcMethod:: r8tx(int NXTLT,int NTHPO,int LENGT,float *CR0,float *CI0)
{
    int J,K;
    float AI0,AI1,AI2,AI3,AI4,AI5,AI6,AI7,
        AR0,AR1,AR2,AR3,AR4,AR5,AR6,AR7,
        BI0,BI1,BI2,BI3,BI4,BI5,BI6,BI7,
        BR0,BR1,BR2,BR3,BR4, BR5,BR6,BR7,
        C1,C2,C3,C4,C5,C6,C7,S1,S2,S3,S4,S5,S6,S7,ARG,SCALE,TI,TR;

    SCALE=PI2/LENGT;
    for(J=0;J<NXTLT;J++)
    {
        ARG=J*SCALE;
        C1=cos(ARG);
        S1=sin(ARG);
        C2=C1*C1-S1*S1;
        S2=C1*S1+C1*S1;
        C3=C1*C2-S1*S2;
        S3=C2*S1+S2*C1;
        C4=C2*C2-S2*S2;
        S4=C2*S2+C2*S2;
        C5=C2*C3-S2*S3;
        S5=C3*S2+S3*C2;
        C6=C3*C3-S3*S3;
        S6=C3*S3+C3*S3;
        C7=C3*C4-S3*S4;
        S7=C4*S3+S4*C3;
        for(K=J;K<NTHPO;K=K+LENGT)
        {
            AR0=CR0[K]+CR0[K+4*NXTLT];
            AR1=CR0[K+NXTLT]+CR0[K+5*NXTLT];
            AR2=CR0[K+2*NXTLT]+CR0[K+6*NXTLT];
            AR3=CR0[K+3*NXTLT]+CR0[K+7*NXTLT];
            AR4=CR0[K]-CR0[K+4*NXTLT];
            AR5=CR0[K+NXTLT]-CR0[K+5*NXTLT];
            AR6=CR0[K+2*NXTLT]-CR0[K+6*NXTLT];
            AR7=CR0[K+3*NXTLT]-CR0[K+7*NXTLT];
            AI0=CI0[K]+CI0[K+4*NXTLT];
            AI1=CI0[K+NXTLT]+CI0[K+5*NXTLT];
            AI2=CI0[K+2*NXTLT]+CI0[K+6*NXTLT];
            AI3=CI0[K+3*NXTLT]+CI0[K+7*NXTLT];
            AI4=CI0[K]-CI0[K+4*NXTLT];
            AI5=CI0[K+NXTLT]-CI0[K+5*NXTLT];
            AI6=CI0[K+2*NXTLT]-CI0[K+6*NXTLT];
            AI7=CI0[K+3*NXTLT]-CI0[K+7*NXTLT];
            BR0=AR0+AR2;
            BR1=AR1+AR3;
            BR2=AR0-AR2;
            BR3=AR1-AR3;
            BR4=AR4-AI6;
            BR5=AR5-AI7;
            BR6=AR4+AI6;
            BR7=AR5+AI7;
            BI0=AI0+AI2;
            BI1=AI1+AI3;
            BI2=AI0-AI2;
            BI3=AI1-AI3;
            BI4=AI4+AR6;
            BI5=AI5+AR7;
            BI6=AI4-AR6;
            BI7=AI5-AR7;
            CR0[K]=BR0+BR1;
            CI0[K]=BI0+BI1;
            if(J > 0)
            {
                CR0[K+NXTLT]=C4*(BR0-BR1)-S4*(BI0-BI1);
                CI0[K+NXTLT]=C4*(BI0-BI1)+S4*(BR0-BR1);
                CR0[K+2*NXTLT]=C2*(BR2-BI3)-S2*(BI2+BR3);
                CI0[K+2*NXTLT]=C2*(BI2+BR3)+S2*(BR2-BI3);
                CR0[K+3*NXTLT]=C6*(BR2+BI3)-S6*(BI2-BR3);
                CI0[K+3*NXTLT]=C6*(BI2-BR3)+S6*(BR2+BI3);
                TR=PI7*(BR5-BI5);
                TI=PI7*(BR5+BI5);
                CR0[K+4*NXTLT]=C1*(BR4+TR)-S1*(BI4+TI);
                CI0[K+4*NXTLT]=C1*(BI4+TI)+S1*(BR4+TR);
                CR0[K+5*NXTLT]=C5*(BR4-TR)-S5*(BI4-TI);
                CI0[K+5*NXTLT]=C5*(BI4-TI)+S5*(BR4-TR);
                TR=-PI7*(BR7+BI7);
                TI=PI7*(BR7-BI7);
                CR0[K+6*NXTLT]=C3*(BR6+TR)-S3*(BI6+TI);
                CI0[K+6*NXTLT]=C3*(BI6+TI)+S3*(BR6+TR);
                CR0[K+7*NXTLT]=C7*(BR6-TR)-S7*(BI6-TI);
                CI0[K+7*NXTLT]=C7*(BI6-TI)+S7*(BR6-TR);
            }
            else
            {
                CR0[K+NXTLT]=BR0-BR1;
                CI0[K+NXTLT]=BI0-BI1;
                CR0[K+2*NXTLT]=BR2-BI3;
                CI0[K+2*NXTLT]=BI2+BR3;
                CR0[K+3*NXTLT]=BR2+BI3;
                CI0[K+3*NXTLT]=BI2-BR3;
                TR=PI7*(BR5-BI5);
                TI=PI7*(BR5+BI5);
                CR0[K+4*NXTLT]=BR4+TR;
                CI0[K+4*NXTLT]=BI4+TI;
                CR0[K+5*NXTLT]=BR4-TR;
                CI0[K+5*NXTLT]=BI4-TI;
                TR=-PI7*(BR7+BI7);
                TI=PI7*(BR7-BI7);
                CR0[K+6*NXTLT]=BR6+TR;
                CI0[K+6*NXTLT]=BI6+TI;
                CR0[K+7*NXTLT]=BR6-TR;
                CI0[K+7*NXTLT]=BI6-TI;
            }
        }
    }


}

void CalcMethod::wavelet( float INTA[1024],float ww[1024],int IP,int ipoint)
{
    /*
    ! 此程序为小波变换
    !*
    !*
    !*        THIS PROGRAM COMPUTE THE ONE-DIMENSIONAL ORTHOGONORMAL
    !*             WAVELETS TRANSFORM OF THE SIGNAL IN vp DATA
    !*               FILE AND RECONSTRUCT THE ORIGINAL SIGNAL
    !*                FROM ITS WAVELET DECOMPOSION USING THE
    !*                     MALLAT'S PYRAMIDAL ALGORITHM
    !*
    !*                           WRITTEN BY liweiyan
    !*                                2004.8
    !*
    !**********************************************************************
    !*
    !*  PARAMETER DISCRIPTION:
    !*
    !*
    !*         N == NUMBER OF DATA SAMPLES TO BE PROCESSED
    !*         M == NUMBER OF fscaleS
    !*         vp(I,J) == THE APPROXIMATION OF THE SIGNAL AT THE
    !*                  RESOLUTION 2^I AND TIME J.(NOTE:vp(0,J)
    !*                  CORRESPONDS TO THE ORIGINAL SIGNAL AT
    !*                  TIME J
    !*         B(I,J) == THE WAVELET DECOMPOSITION RESULT AT
    !*                   RESOLUTION 2^I AND TIME J.(NOTE:I>=1)
    !*         AA(I,J) == THE RECONSTRUCTED SIGNAL AT THE RESOLUTION
    !*                   2^I AND TIME J.(NOTE:I>=0)
    !*         fscale(I) == THE WAVELET COMPOSITION AT RESPONDING
    !*                     fscale IS VALUED 0
    !*
    !*  SUBROUTINE:
    !*
    !*         DWT ==  COMPUTE THE ONE-DIMENSIONAL ORTOGONORMAL
    !*                 WAVELET TRANSFORM
    !*         RECON_DWT == RECONSTRUCT THE ORIGINAL SIGNAL FROM
    !*                      ITS WAVELET DECOMPOSITION
    !*
    !*****************************************************************
    */


    //!   输入参数名，个数与moC对应
    //   ! 下面可以定义处理程序的变量
    float vp[7][1024],B[7][1024],AA[7][1024];
    float C1[8],H[4],G[2],K1[6];
    float wmax,wwmax;
    int M,N,fscale[7],i,k,I512;

    C1[0]=1.5;
    C1[1]=1.12;
    C1[2]=1.03;
    C1[3]=1.01;
    C1[4]=1.;
    C1[5]=1.;
    C1[6]=1.;
    C1[7]=1.;
    H[0]=0.125;
    H[1]=0.375;
    H[2]=0.375;
    H[3]=0.125;
    G[0]=-2.0;
    G[1]=2.0;

    K1[0]=0.0078125;
    K1[1]=0.054685;
    K1[2]=0.171875;
    K1[3]=-0.171875;
    K1[4]=-0.054685;
    K1[5]=-0.0078125;


    M=5;
    fscale[1]=0;
    fscale[2]=0;
    fscale[3]=0;
    fscale[4]=0;
    fscale[0]=0;

    if(IP<0)
    {
        fscale[1]=1;
        fscale[2]=1;
        fscale[3]=1;
        fscale[4]=0;
        fscale[0]=1;
        //	fscale[IP]=1;
    }
    else
    {
        fscale[IP-1]=1;
        fscale[IP]=1;
    }
    I512=ipoint;
    N=512;
    if(I512>512)N=1024;

    wmax= vp[0][0]=fabs(INTA[0]);


    for(i=0;i<ipoint;i++)
    {
        vp[0][i]=INTA[i];
        if(fabs(vp[0][i])>wmax)wmax=fabs(vp[0][i]);
    }
    if(ipoint<N)
    {
        for( i=ipoint;i<N;i++)vp[0][i]=0.0;

    }
    DWT(vp,B,M,N,H,C1,G);
    for(i=0;i<N+10;i++)AA[M][i]=vp[M][i];

    for( k=0;k<M;k++)
    {
        if(fscale[k]==0)
        {
            for( i=0;i<N+10;i++)B[k][i]=0.0;
        }
    }
    RECON_DWT(AA,B,M,N,H,C1,K1);
    wwmax=fabs(AA[0][0]);
    for(i=0;i<ipoint;i++)
    {
        vp[0][i]=AA[0][i];
        if(fabs(vp[0][i])>wwmax)wwmax=fabs(vp[0][i]);
    }
    for( i=0;i<N-1;i++)ww[i]=AA[0][i]*wmax/wwmax;
    return;
}

/*
!**
!**      COMPUTE THE ONE-DIMENSIONAL ORTHOGONAOTMAL
!**       WAVELETS TRANSFORM OF THE SIGNAL IN THE
!**                   DATA FILE FWAVE
!**
!**
*/
void CalcMethod::DWT(float vp[7][1024],float B[7][1024],int M,int N,float H[4],float C1[8],float G[2])
{

    int I,J,K,C;
    float P,Q,C2;
    C=1;
    for( K=0;K<M;K++)
    {
        C2=1.0/C1[K-1];
        if(K==0)
        {
            for( I=N;I<N+2*C+1;I++)vp[K][I]=vp[K][2*N-1-I];

            for( I=0;I<2;I++)
            {
                P=0.0;
                Q=0.0;
                for( J=-1;J<=2;J++)P=P+H[J+1]*vp[0][int(fabs(I-J+0.5)-0.5)];
                for( J=0;J<=1;J++)Q=Q+G[J]*vp[0][int(fabs(I-J+0.5)-0.5)];

                vp[0][I]=P;
                B[0][I]=-C2*Q;
            }

            for( I=3;I<N;I++)
            {
                P=0.0;
                Q=0.0;
                for( J=-1;J<=2;J++)P=P+H[J+1]*vp[0][I-J];
                for( J=0;J<=1;J++)Q=Q+G[J]*vp[0][I-J];

                vp[0][I]=P;
                B[0][I]=-C2*Q;
            }

        }
        else
        {
            for(I=N+1;I<=N+2*C+1;I++) vp[K-1][I]=vp[K-1][2*N-I];

            for(I=C/2;I<=2*C;I++)
            {
                P=0.0;
                Q=0.0;
                for(J=-1;J<=2;J++)P=P+H[J+1]*vp[K-1][int(fabs((double)(I-C*J)))];

                for(J=0;J<=1;J++)  Q=Q+G[J]*vp[K-1][int(fabs((double)(I-C*J)))];

                vp[K][I-C/2]=P;
                B[K][I-C/2]=-C2*Q;
            }
            for( I=2*C+1;I<=N+C/2;I++)
            {
                P=0.0;
                Q=0.0;
                for(J=-1;J<=2;J++)P=P+H[J+1]*vp[K-1][I-C*J];

                for(J=0;J<=1;J++)  Q=Q+G[J]*vp[K-1][I-C*J];

                vp[K][I-C/2]=P;
                B[K][I-C/2]=-C2*Q;
            }
        }
        C=2*C;
    }
    return;
}
/*
!**
!**       RECONSTRUCT THE ORIGINAL SIGNAL FROM ITS
!**         WAVELET DECOMPOSITION B
!**
*/
void CalcMethod:: RECON_DWT(float AA[7][1024],float B[7][1024], int M,int N,float H[4],float C1[8],float K1[6])

{
    int I,J,K,C;
    float P,Q,C2;
    C=pow(2.0,(M-1));
    for( K=M-1;K>0;K=K-1)
    {
        C2=C1[K-1];
        for( I=N+1;I<=N+3*C;I++)
        {
            B[K][I]=-B[K][2*N-I];
            AA[K][I]=AA[K][2*N-I];
        }
        for( I=-C/2;I<=3*C;I++)
        {
            P=0.0;
            Q=0.0;
            for( J=-1;J<=2;J++) P=P+H[J+1]*AA[K][int(fabs((double)(I+C*J)))];

            for( J=-3;J<=2;J++)
            {
                if((I-C*J)<0)
                {
                    Q=Q+K1[J+3]*(-1)*B[K][int(fabs((double)(I-C*J)))];
                }
                else
                {
                    Q=Q+K1[J+3]*B[K][I-C*J];
                }
            }
            AA[K-1][I+C/2]=P-C2*Q;
        }
        for(I=3*C+1;I<=N-C/2;I++)
        {
            P=0.0;
            Q=0.0;
            for( J=-1;J<=2;J++) P=P+H[J+1]*AA[K][I+C*J];

            for( J=-3;J<=2;J++) Q=Q+K1[J+3]*B[K][I-C*J];

            AA[K-1][I+C/2]=P-C2*Q;
        }
        C=C/2;
    }
    return;
}


complex<float>  CalcMethod::GetRootMonteCarlo(complex<float> k0,float vp,float vs,float vf,float denf,float den,float R,float w,int n,float MTOOL,float R_TOOL,float df)
{
    int mm=0;
    int k;
    float xx,yy,a,r,z,x1,y1,z1,eps=1.0e-10,k_real,k_imag;
    int nControlB=10;
    float xStart;
    //	float PI=3.1415926;
    complex<float> kw0,kw;

    //	w=2*PI*fff;
    //	k0=complex<float>(w/vf/0.85,0.0);
    if(n==0)
    {
        if(vs>vf)
        {
            kw0=k0+complex <float>(0.1,0.);//complex <float>(w/vf+.10,0.0);
        }
        else
        {
            if(w>2000*2*3.1415926)
            {
                kw0=k0-complex <float>(0.1,0.);//complex <float>(w/vf-.10,0.0);
            }
            else
            {
                kw0=k0+complex <float>(0.1,0.);
            }
        }



        //	 kw0=initialvalue(real(k0),imag(k0),vp,vs,vf,denf,den,R,w,n,MTOOL,R_TOOL,df);
        xStart=0.1;
    }//w/vf;
    if(n==1)
    {
        kw0=complex<float>(w/vs/0.85,0.0);
        //	 kw0=initialvalue(real(k0),imag(k0),vp,vs,vf,denf,den,R,w,n,MTOOL,R_TOOL,df);
        xStart=.10;//w/vf;
    }
    a=xStart;
    k=1;
    r=1.0;
    xx=real(kw0);
    yy=imag(kw0);
    z=Func(xx,yy,vp,vs,vf,denf,den,R, w,n,MTOOL,R_TOOL);// whp i=3时，kew[i]=CalcMethod::GetRootMonteCarlo(k0,

    // ???è?????ó?vp
    while (a>=eps)
    {
        x1=-a+2.0*a*rnd(&r);
        x1=xx+x1;
        y1=-a+2.0*a*rnd(&r);
        y1=yy+y1;

        z1=Func(x1,y1,vp,vs,vf,denf,den,R,w,n,MTOOL,R_TOOL);

        k=k+1;
        if (z1>=z)
        {
            if (k>nControlB)
            {
                k=1;
                a=a/2.0;
                mm++;
            }
        }
        else
        {
            //			k=1;
            xx=x1;
            yy=y1;
            z=z1;
            if (z<eps)
            {
                k_real=xx;
                k_imag=yy;
                kw=complex <float>(xx,yy);
                return kw;
            }
        }
    }

    k_real=xx;
    k_imag=yy;
    kw=complex <float>(xx,yy);
    return kw;
}




float CalcMethod::rnd(float* r)
{
    int m;
    float s,u,v,p;

    s=65536.0;
    u=2053.0;
    v=13849.0;
    m=(int)(*r/s);
    *r=*r-m*s;
    *r=u*(*r)+v;
    m=(int)(*r/s);
    *r=*r-m*s; p=*r/s;

    return(p);
}



float CalcMethod::Func(float x,float y ,float vp,float vs,float vf,float denf,float den,float R,float w,int n,float MTOOL,float R_TOOL)
{
    complex<float>  M[4][4],C,fun,p,s,f,g,k,ttt;
    complex<float>fr,sr,pr,fa,ai,inter[4],Inter;
    float arg0=0.0,SPACE=3.6;
    float kp,kf,ks,etool,temp;
    int i,j;
    int u[24][4]={{0,1,2,3},{0,1,3,2},{0,2,1,3},{0,2,3,1},{0,3,1,2},{0,3,2,1},
    {1,0,2,3},{1,0,3,2},{1,2,0,3},{1,2,3,0},{1,3,0,2},{1,3,2,0},
    {2,0,3,1},{2,0,1,3},{2,1,3,0},{2,1,0,3},{2,3,0,1},{2,3,1,0},
    {3,0,1,2},{3,0,2,1},{3,1,2,0},{3,1,0,2},{3,2,0,1},{3,2,1,0}};
    int t,sd[4];
    //       k=w/v;
    //		w=w/1000.;
    //		ttt=complex<float>(1.e-10,0.0);
    k=complex<float>(x,y);
    kp=w/vp;
    ks=w/vs;
    kf=w/vf;
    ai=complex<float>(0.0,1.0);
    ttt=(k-complex<float>(kf,0.0))*(k-complex<float>(kf,0.0))*R*R;
    p=sqrt(k*k-complex<float>(kp*kp,0));
    s=sqrt(k*k-complex<float>(ks*ks,0));
    f=sqrt(k*k-complex<float>(kf*kf,0));
    complex<float> common;
    inter[0]=pow(f*complex<float>(R_TOOL/2,0),n);
    inter[1]=complex<float>(cos(n*(arg0-3.1415926/2)),0);
    Inter=ai*k*complex<float>(SPACE,0);
    inter[2]=complex<float>(pow(2.718282,(double)real(Inter))*cos(imag(Inter)),pow(2.718282,(double)real(Inter))*sin(imag(Inter)));
    common=inter[3]=complex<float>(1.0,0.0);//inter[0]*inter[1]*inter[2]/complex<float>(1000.,0.0);


    fr=f*complex <float>(R,0.0);
    pr=p*complex <float>(R,0.0);
    sr=s*complex <float>(R,0.0);
    fa=f*complex <float>(R_TOOL,0.0);

    temp=MTOOL/R_TOOL;
    etool=(temp*real(f*Ibessel(1,fa))+denf*w*w*real(Ibessel(0,fa)))/(temp*real(f*Kbessel(1,fa))-denf*w*w*real(Kbessel(0,fa)));
    M[0][0]=complex <float>(-1.,0.0)*complex <float>(n/R,0.0)*Ibessel(n,fr)-f*Ibessel(n+1,fr);
    M[0][1]=complex <float>(-1.,0.0)*p*(complex <float>(-n,0)/(pr)*Kbessel(n,pr)+Kbessel(n+1,pr));
    M[0][2]=complex <float>(-1.,0.0)*complex <float>(n/R,0)*Kbessel(n,sr);
    M[0][3]=ai*complex <float>(-ks,0.0)*(complex <float>(-n,0.0)/(sr)*Kbessel(n,sr)+Kbessel(n+1,sr));
    M[1][0]=(complex <float>(denf*w*w,0.0)*Ibessel(n,fr));
    M[1][1]=(complex <float>(den,0)*(complex <float>(2.0*vs*vs,0)*k*k-complex <float>(w*w,0.0))*Kbessel(n,pr)+complex <float>(2.0*den*vs*vs/R,0.0)*p*((complex <float>(n*(n-1),0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[1][2]=(complex <float>(-2*n*den*vs*vs/R,0.0)*s*((complex <float>(1-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    M[1][3]=(ai*k*complex <float>(2*den*vs*vs,0.0)*s*s*Kbessel(n,sr)+ai*complex <float>(2*ks*den*vs*vs/R,0)*((complex <float>(n*(n-1),0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    M[2][0]=(0.0,0.0);
    M[2][1]=(complex <float>(-2*n*den*vs*vs,0.0)*p*((complex <float>(1-n,0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[2][2]=(complex <float>(-den*vs*vs,0)*s*s*((complex <float>(1.0,0.0)+complex <float>(2*n*(n-1),0.0)/(sr*sr))*Kbessel(n,sr)+(complex <float>(2,0.0)/sr)*Kbessel(n+1,sr)));
    M[2][3]=(ai*complex <float>(2*ks*n*den*vs*vs/R,0)*((complex <float>(1-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    M[3][0]=(0.0,0.0);
    M[3][1]=(ai*k*complex <float>(-2*den*vs*vs,0)*p*((complex <float>(-n,0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[3][2]=(ai*k*complex <float>(-n*den*vs*vs/R,0)*Kbessel(n,sr));
    M[3][3]=((k*k+s*s)*s*complex <float>(den*vs*vs,0)*((complex <float>(-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    if(MTOOL>0.01)
    {
        if(n==0)
        {
            complex<float> aaa=Ibessel(n,fr);
            aaa=Ibessel(n+1,fr);
            aaa=Kbessel(n+1,fr);//whp 溢出
            M[0][0]=complex <float>(-1.,0.0)*complex <float>(n/R,0.0)*Ibessel(n,fr)-f*(Ibessel(n+1,fr)-etool*Kbessel(n+1,fr));
            M[1][0]=(complex <float>(denf*w*w,0.0)*(Ibessel(n,fr)+etool*Kbessel(n,fr)));
        }
        if(n==1)
        {
            M[0][0]=complex <float>(-1.,0.0)*complex <float>(n/R,0.0)*(Ibessel(n,fr)-etool*Kbessel(n,fr))-f*(Ibessel(n+1,fr));
            M[1][0]=(complex <float>(denf*w*w,0.0)*(Ibessel(n,fr)-etool*Kbessel(n,fr)));
        }

    }
    for(j=1;j<=3;j++)
    {
        for(i=0;i<=3;i++)M[j][i]=M[j][i]*ttt;
    }

    for(j=0;j<24;j++)
    {
        sd[0]=u[j][0];
        sd[1]=u[j][1];
        sd[2]=u[j][2];
        sd[3]=u[j][3];

        t=reversion(sd);
        //			complex <float> temp1=complex <float>(pow(-1,t),0.0);
        //            complex <float> temp2=complex <float>((-1)^t,0.0);
        fun=fun+complex <float>(pow(-1.0,t),0.0)*M[0][u[j][0]]*common*M[1][u[j][1]]*common*M[2][u[j][2]]*common*M[3][u[j][3]]*common;

    }
    //whp add 2010.9.15 for sqrt(real(fun)*real(fun)+imag(fun)*imag(fun));出现浮点溢出
    float aa=fabs(real(fun));
    float bb=fabs(imag(fun));
    float mud=1;
    while(aa>=1e+18||bb>=1e+18)
    {
        aa/=10000.;bb/=10000.;mud*=100.;
    }
    aa=aa*aa;
    bb=bb*bb;
    float cc=aa+bb;
    cc=sqrt(cc)*mud;
    //float aaa=sqrt(real(fun)*real(fun)+imag(fun)*imag(fun));
    return cc;//sqrt(real(fun)*real(fun)+imag(fun)*imag(fun));
}

//
float CalcMethod::FuncD(float x,float y ,float vp,float vs,float vf,float denf,float den,float R,float w,int n,float MTOOL,float R_TOOL)
{
    complex<float>  M[4][4],C,fun,p,s,f,g,k,ttt;
    complex<float>fr,sr,pr,fa,ai,inter[4],Inter;
    float arg0=0.0,SPACE=3.6,Funcd;
    float kp,kf,ks,etool,temp,E10=1000000000000.;
    int i,j;
    int u[24][4]={{0,1,2,3},{0,1,3,2},{0,2,1,3},{0,2,3,1},{0,3,1,2},{0,3,2,1},
    {1,0,2,3},{1,0,3,2},{1,2,0,3},{1,2,3,0},{1,3,0,2},{1,3,2,0},
    {2,0,3,1},{2,0,1,3},{2,1,3,0},{2,1,0,3},{2,3,0,1},{2,3,1,0},
    {3,0,1,2},{3,0,2,1},{3,1,2,0},{3,1,0,2},{3,2,0,1},{3,2,1,0}};
    int t,sd[4];
    //       k=w/v;
    //		w=w/1000.;
    //		ttt=complex<float>(1.e-10,0.0);
    k=complex<float>(x,y);
    kp=w/vp;
    ks=w/vs;
    kf=w/vf;
    ai=complex<float>(0.0,1.0);
    ttt=(k-complex<float>(kf,0.0))*(k-complex<float>(kf,0.0))*R*R;
    p=sqrt(k*k-complex<float>(kp*kp,0));
    s=sqrt(k*k-complex<float>(ks*ks,0));
    f=sqrt(k*k-complex<float>(kf*kf,0));
    complex<float> common;
    inter[0]=pow(f*complex<float>(R_TOOL/2,0),n);
    inter[1]=complex<float>(cos(n*(arg0-3.1415926/2)),0);
    Inter=ai*k*complex<float>(SPACE,0);
    inter[2]=complex<float>(pow(2.718282,(double)real(Inter))*cos(imag(Inter)),pow(2.718282,(double)real(Inter))*sin(imag(Inter)));
    common=inter[3]=complex<float>(1.0,0.0);//inter[0]*inter[1]*inter[2]/complex<float>(1000.,0.0);


    fr=f*complex <float>(R,0.0);
    pr=p*complex <float>(R,0.0);
    sr=s*complex <float>(R,0.0);
    fa=f*complex <float>(R_TOOL,0.0);

    temp=MTOOL/R_TOOL;
    etool=(temp*real(f*Ibessel(1,fa))+denf*w*w*real(Ibessel(0,fa)))/(temp*real(f*Kbessel(1,fa))-denf*w*w*real(Kbessel(0,fa)));
    M[0][0]=complex <float>(-1.,0.0)*complex <float>(n/R,0.0)*Ibessel(n,fr)-f*Ibessel(n+1,fr);
    M[0][1]=complex <float>(-1.,0.0)*p*(complex <float>(-n,0)/(pr)*Kbessel(n,pr)+Kbessel(n+1,pr));
    M[0][2]=complex <float>(-1.,0.0)*complex <float>(n/R,0)*Kbessel(n,sr);
    M[0][3]=ai*complex <float>(-ks,0.0)*(complex <float>(-n,0.0)/(sr)*Kbessel(n,sr)+Kbessel(n+1,sr));
    M[1][0]=(complex <float>(denf*w*w,0.0)*Ibessel(n,fr));
    M[1][1]=(complex <float>(den,0)*(complex <float>(2.0*vs*vs,0)*k*k-complex <float>(w*w,0.0))*Kbessel(n,pr)+complex <float>(2.0*den*vs*vs/R,0.0)*p*((complex <float>(n*(n-1),0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[1][2]=(complex <float>(-2*n*den*vs*vs/R,0.0)*s*((complex <float>(1-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    M[1][3]=(ai*k*complex <float>(2*den*vs*vs,0.0)*s*s*Kbessel(n,sr)+ai*complex <float>(2*ks*den*vs*vs/R,0)*((complex <float>(n*(n-1),0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    M[2][0]=(0.0,0.0);
    M[2][1]=(complex <float>(-2*n*den*vs*vs,0.0)*p*((complex <float>(1-n,0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[2][2]=(complex <float>(-den*vs*vs,0)*s*s*((complex <float>(1.0,0.0)+complex <float>(2*n*(n-1),0.0)/(sr*sr))*Kbessel(n,sr)+(complex <float>(2,0.0)/sr)*Kbessel(n+1,sr)));
    M[2][3]=(ai*complex <float>(2*ks*n*den*vs*vs/R,0)*((complex <float>(1-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    M[3][0]=(0.0,0.0);
    M[3][1]=(ai*k*complex <float>(-2*den*vs*vs,0)*p*((complex <float>(-n,0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[3][2]=(ai*k*complex <float>(-n*den*vs*vs/R,0)*Kbessel(n,sr));
    M[3][3]=((k*k+s*s)*s*complex <float>(den*vs*vs,0)*((complex <float>(-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    if(MTOOL>0.01)
    {
        if(n==0)
        {
            M[0][0]=complex <float>(-1.,0.0)*complex <float>(n/R,0.0)*Ibessel(n,fr)-f*(Ibessel(n+1,fr)-etool*Kbessel(n+1,fr));
            M[1][0]=(complex <float>(denf*w*w,0.0)*(Ibessel(n,fr)+etool*Kbessel(n,fr)));
        }
        if(n==1)
        {
            M[0][0]=complex <float>(-1.,0.0)*complex <float>(n/R,0.0)*(Ibessel(n,fr)-etool*Kbessel(n,fr))-f*(Ibessel(n+1,fr));
            M[1][0]=(complex <float>(denf*w*w,0.0)*(Ibessel(n,fr)-etool*Kbessel(n,fr)));
        }

    }
    for(j=1;j<=3;j++)
    {
        for(i=0;i<=3;i++)M[j][i]=M[j][i]*ttt;
    }

    for(j=0;j<24;j++)
    {
        sd[0]=u[j][0];
        sd[1]=u[j][1];
        sd[2]=u[j][2];
        sd[3]=u[j][3];

        t=reversion(sd);
        //			complex <float> temp1=complex <float>(pow(-1,t),0.0);
        //            complex <float> temp2=complex <float>((-1)^t,0.0);
        fun=fun+complex <float>(pow(-1.0,t),0.0)*M[0][u[j][0]]*common*M[1][u[j][1]]*common*M[2][u[j][2]]*common*M[3][u[j][3]]*common;

    }
    Funcd= real(fun)/E10;
    return Funcd;
}
//
float CalcMethod::Func_M11(float x,float y ,float vp,float vs,float vf,float denf,float den,float R,float w,int n,float MTOOL,float R_TOOL)
{
    complex<float>  M[3][3],C,fun,p,s,f,g,k,ttt;
    complex<float>fr,sr,pr,fa,ai,inter[4],Inter;
    float arg0=0.0,SPACE=3.6;
    float kp,kf,ks,etool,temp,m11,E10=1000000000000.;
    int i,j;

    int u[6][3]={{0,1,2},{0,2,1},
    {1,0,2},{1,2,0},
    {2,0,1},{2,1,0}};
    int t,sd[3];
    //       k=w/v;
    //		w=w/1000.;
    //		ttt=complex<float>(1.e-10,0.0);
    k=complex<float>(x,y);
    kp=w/vp;
    ks=w/vs;
    kf=w/vf;
    ai=complex<float>(0.0,1.0);
    ttt=(k-complex<float>(kf,0.0))*(k-complex<float>(kf,0.0))*R*R;
    p=sqrt(k*k-complex<float>(kp*kp,0));
    s=sqrt(k*k-complex<float>(ks*ks,0));
    f=sqrt(k*k-complex<float>(kf*kf,0));
    complex<float> common;
    inter[0]=pow(f*complex<float>(R_TOOL/2,0),n);
    inter[1]=complex<float>(cos(n*(arg0-3.1415926/2)),0);
    Inter=ai*k*complex<float>(SPACE,0);
    inter[2]=complex<float>(pow(2.718282,(double)real(Inter))*cos(imag(Inter)),pow(2.718282,(double)real(Inter))*sin(imag(Inter)));
    common=inter[3]=complex<float>(1.0,0.0);//inter[0]*inter[1]*inter[2]/complex<float>(1000.,0.0);


    fr=f*complex <float>(R,0.0);
    pr=p*complex <float>(R,0.0);
    sr=s*complex <float>(R,0.0);
    fa=f*complex <float>(R_TOOL,0.0);

    temp=MTOOL/R_TOOL;
    etool=(temp*real(f*Ibessel(1,fa))+denf*w*w*real(Ibessel(0,fa)))/(temp*real(f*Kbessel(1,fa))-denf*w*w*real(Kbessel(0,fa)));
    //		M[0][0]=complex <float>(-1.,0.0)*complex <float>(n/R,0.0)*Ibessel(n,fr)-f*Ibessel(n+1,fr);
    //		M[0][1]=complex <float>(-1.,0.0)*p*(complex <float>(-n,0)/(pr)*Kbessel(n,pr)+Kbessel(n+1,pr));
    //		M[0][2]=complex <float>(-1.,0.0)*complex <float>(n/R,0)*Kbessel(n,sr);
    //		M[0][3]=ai*complex <float>(-ks,0.0)*(complex <float>(-n,0.0)/(sr)*Kbessel(n,sr)+Kbessel(n+1,sr));
    //		M[1][0]=(complex <float>(denf*w*w,0.0)*Ibessel(n,fr));
    M[0][0]=(complex <float>(den,0)*(complex <float>(2.0*vs*vs,0)*k*k-complex <float>(w*w,0.0))*Kbessel(n,pr)+complex <float>(2.0*den*vs*vs/R,0.0)*p*((complex <float>(n*(n-1),0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[0][1]=(complex <float>(-2*n*den*vs*vs/R,0.0)*s*((complex <float>(1-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    M[0][2]=(ai*k*complex <float>(2*den*vs*vs,0.0)*s*s*Kbessel(n,sr)+ai*complex <float>(2*ks*den*vs*vs/R,0)*((complex <float>(n*(n-1),0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    //		M[2][0]=(0.0,0.0);
    M[1][0]=(complex <float>(-2*n*den*vs*vs,0.0)*p*((complex <float>(1-n,0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[1][1]=(complex <float>(-den*vs*vs,0)*s*s*((complex <float>(1.0,0.0)+complex <float>(2*n*(n-1),0.0)/(sr*sr))*Kbessel(n,sr)+(complex <float>(2,0.0)/sr)*Kbessel(n+1,sr)));
    M[1][2]=(ai*complex <float>(2*ks*n*den*vs*vs/R,0)*((complex <float>(1-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    //		M[3][0]=(0.0,0.0);
    M[2][0]=(ai*k*complex <float>(-2*den*vs*vs,0)*p*((complex <float>(-n,0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[2][1]=(ai*k*complex <float>(-n*den*vs*vs/R,0)*Kbessel(n,sr));
    M[2][2]=((k*k+s*s)*s*complex <float>(den*vs*vs,0)*((complex <float>(-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));

    for(j=0;j<=2;j++)
    {
        for(i=0;i<=2;i++)M[j][i]=M[j][i]*ttt;
    }

    for(j=0;j<6;j++)
    {
        sd[0]=u[j][0];
        sd[1]=u[j][1];
        sd[2]=u[j][2];
        //			sd[3]=u[j][3];

        t=reversion3(sd);

        fun=fun+complex <float>(pow(-1.0,t),0.0)*M[0][u[j][0]]*common*M[1][u[j][1]]*common*M[2][u[j][2]]*common;

    }
    m11=real(fun)/E10;
    return m11;
}
//
//
float CalcMethod::Func_M21(float x,float y ,float vp,float vs,float vf,float denf,float den,float R,float w,int n,float MTOOL,float R_TOOL)
{
    complex<float>  M[3][3],C,fun,p,s,f,g,k,ttt;
    complex<float>fr,sr,pr,fa,ai,inter[4],Inter;
    float arg0=0.0,SPACE=3.6,m21;
    float kp,kf,ks,etool,temp,E10=1000000000000.;
    int i,j;
    int u[6][3]={{0,1,2},{0,2,1},
    {1,0,2},{1,2,0},
    {2,0,1},{2,1,0}};
    int t,sd[3];
    //       k=w/v;
    //		w=w/1000.;
    //		ttt=complex<float>(1.e-10,0.0);
    k=complex<float>(x,y);
    kp=w/vp;
    ks=w/vs;
    kf=w/vf;
    ai=complex<float>(0.0,1.0);
    ttt=(k-complex<float>(kf,0.0))*(k-complex<float>(kf,0.0))*R*R;
    p=sqrt(k*k-complex<float>(kp*kp,0));
    s=sqrt(k*k-complex<float>(ks*ks,0));
    f=sqrt(k*k-complex<float>(kf*kf,0));
    complex<float> common;
    inter[0]=pow(f*complex<float>(R_TOOL/2,0),n);
    inter[1]=complex<float>(cos(n*(arg0-3.1415926/2)),0);
    Inter=ai*k*complex<float>(SPACE,0);
    inter[2]=complex<float>(pow(2.718282,(double)real(Inter))*cos(imag(Inter)),pow(2.718282,(double)real(Inter))*sin(imag(Inter)));
    common=inter[3]=complex<float>(1.0,0.0);//inter[0]*inter[1]*inter[2]/complex<float>(1000.,0.0);


    fr=f*complex <float>(R,0.0);
    pr=p*complex <float>(R,0.0);
    sr=s*complex <float>(R,0.0);
    fa=f*complex <float>(R_TOOL,0.0);

    temp=MTOOL/R_TOOL;
    etool=(temp*real(f*Ibessel(1,fa))+denf*w*w*real(Ibessel(0,fa)))/(temp*real(f*Kbessel(1,fa))-denf*w*w*real(Kbessel(0,fa)));
    //		M[0][0]=complex <float>(-1.,0.0)*complex <float>(n/R,0.0)*Ibessel(n,fr)-f*Ibessel(n+1,fr);
    M[0][0]=complex <float>(-1.,0.0)*p*(complex <float>(-n,0)/(pr)*Kbessel(n,pr)+Kbessel(n+1,pr));
    M[0][1]=complex <float>(-1.,0.0)*complex <float>(n/R,0)*Kbessel(n,sr);
    M[0][2]=ai*complex <float>(-ks,0.0)*(complex <float>(-n,0.0)/(sr)*Kbessel(n,sr)+Kbessel(n+1,sr));
    //		M[1][0]=(complex <float>(denf*w*w,0.0)*Ibessel(n,fr));
    //        M[1][1]=(complex <float>(den,0)*(complex <float>(2.0*vs*vs,0)*k*k-complex <float>(w*w,0.0))*Kbessel(n,pr)+complex <float>(2.0*den*vs*vs/R,0.0)*p*((complex <float>(n*(n-1),0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    //		M[1][2]=(complex <float>(-2*n*den*vs*vs/R,0.0)*s*((complex <float>(1-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    //		M[1][3]=(ai*k*complex <float>(2*den*vs*vs,0.0)*s*s*Kbessel(n,sr)+ai*complex <float>(2*ks*den*vs*vs/R,0)*((complex <float>(n*(n-1),0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    //		M[2][0]=(0.0,0.0);
    M[1][0]=(complex <float>(-2*n*den*vs*vs,0.0)*p*((complex <float>(1-n,0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[1][1]=(complex <float>(-den*vs*vs,0)*s*s*((complex <float>(1.0,0.0)+complex <float>(2*n*(n-1),0.0)/(sr*sr))*Kbessel(n,sr)+(complex <float>(2,0.0)/sr)*Kbessel(n+1,sr)));
    M[1][2]=(ai*complex <float>(2*ks*n*den*vs*vs/R,0)*((complex <float>(1-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    //		M[3][0]=(0.0,0.0);
    M[2][0]=(ai*k*complex <float>(-2*den*vs*vs,0)*p*((complex <float>(-n,0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[2][1]=(ai*k*complex <float>(-n*den*vs*vs/R,0)*Kbessel(n,sr));
    M[2][2]=((k*k+s*s)*s*complex <float>(den*vs*vs,0)*((complex <float>(-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));

    for(j=0;j<=2;j++)
    {
        for(i=0;i<=2;i++)M[j][i]=M[j][i]*ttt;
    }

    for(j=0;j<6;j++)
    {
        sd[0]=u[j][0];
        sd[1]=u[j][1];
        sd[2]=u[j][2];
        //			sd[3]=u[j][3];

        t=reversion3(sd);

        fun=fun+complex <float>(pow(-1.0,t),0.0)*M[0][u[j][0]]*common*M[1][u[j][1]]*common*M[2][u[j][2]]*common;

    }
    m21=real(fun)/E10;
    return m21;
}
//
float CalcMethod::Get_EW(float k0,float vp,float vs,float vf,float denf,float den,float R,float w,int n,float MTOOL,float R_TOOL,float df)
{
    complex<float> uff,pff,p,s,f,g,k,ttt;
    complex<float>fr,sr,pr,fa,ai;
    float M11,M21,arg0=0.0,SPACE=3.6;
    float kp,kf,ks,x,dx;

    ai=(0.,1.);
    //       k=w/v;
    //		w=w/1000.;
    //		ttt=complex<float>(1.e-10,0.0);
    k=complex<float>(k0,0);
    kp=w/vp;
    ks=w/vs;
    kf=w/vf;
    ai=complex<float>(0.0,1.0);
    ttt=(k-complex<float>(kf,0.0))*(k-complex<float>(kf,0.0))*R*R;
    p=sqrt(k*k-complex<float>(kp*kp,0));
    s=sqrt(k*k-complex<float>(ks*ks,0));
    f=sqrt(k*k-complex<float>(kf*kf,0));


    fr=f*complex <float>(R,0.0);
    pr=p*complex <float>(R,0.0);
    sr=s*complex <float>(R,0.0);
    fa=f*complex <float>(R_TOOL,0.0);
    uff=-f*Ibessel(1,fr)*exp(ai*k*SPACE);
    pff=-denf*w*w*Ibessel(1,fr)*exp(ai*k*SPACE);
    M11=Func_M11(k0,0 , vp, vs, vf, denf, den, R, w,n, MTOOL, R_TOOL);
    M21=Func_M21(k0,0 , vp, vs, vf, denf, den, R, w,n, MTOOL, R_TOOL);
    dx=FuncD(k0+0.001,0 , vp, vs, vf, denf, den, R, w,n, MTOOL, R_TOOL)-FuncD(k0,0 , vp, vs, vf, denf, den, R, w,n, MTOOL, R_TOOL);
    x=dx/0.001;

    return real(uff*M11-pff*M21)/x;
}
//
complex<float> CalcMethod::Ibessel(int n,complex<float> x)
{
    complex<float> k,wi[2],wk[2];
    int IONLY=0;
    ///	if(fabs(imag(x))<=0.0000000001)
    ///	{
    ///        return complex <float>(TransformativeIntegerBessel1stFunction(n,real(x)),0.0);
    ///	}
    ///	else
    {
        if(n==0)
        {
            CFB(x, wi, wk,IONLY);
            return wi[0];
        }
        if(n==1)
        {
            CFB(x, wi, wk, IONLY);
            return wi[1];
        }
        k=(complex<float>(-2*n,0)/x)*Ibessel(n-1,x)+Ibessel(n-2,x);
        return k;
    }
}

complex<float> CalcMethod::Kbessel(int n,complex<float> x)
{
    complex<float> k,wi[2],wk[2];
    int IONLY=0;
    ///    if(fabs(imag(x))<=0.0000000001)
    ///	{
    ///       return complex <float>(TransformativeIntegerBessel2ndFunction(n,real(x)),0.0);
    ///	}
    ///	else
    ///	{
    if(n==0)
    {
        CFB(x, wi, wk, IONLY);

        return wk[0];
    }
    if(n==1)
    {
        CFB(x, wi, wk, IONLY);
        return wk[1];
    }
    k=(complex<float>(2*n,0.0)/x)*Ibessel(n-1,x)+Ibessel(n-2,x);
    return k;
    ///	}
}
inline bool   CalcMethod::doubleEqual(float lhs, float rhs)
{
    float	DOUBLEERROR=0.000001;
    if (fabs(lhs - rhs) < DOUBLEERROR)
        return true;
    else
        return false;
}

int CalcMethod::reversion(int vp[4])
{
    int i,t=0,j;
    for(i=1;i<4;i++)
    {
        for(j=0;j<i;j++)
        {
            if(vp[i]<vp[j])
            {
                t=t+1;
            }
        }

    }
    return t;
}
int CalcMethod::reversion3(int vp[3])
{
    int i,t=0,j;
    for(i=1;i<3;i++)
    {
        for(j=0;j<i;j++)
        {
            if(vp[i]<vp[j])
            {
                t=t+1;
            }
        }
    }
    return t;
}
float CalcMethod::TransformativeIntegerBessel1stFunction(int n,float x)
{
    int i,m;
    float t,y,p,b0,b1,q;
    float a[7]={ 1.0,3.5156229,3.0899424,1.2067492,
        0.2659732,0.0360768,0.0045813};
    float b[7]={ 0.5,0.87890594,0.51498869,
        0.15084934,0.02658773,0.00301532,0.00032411};
    float c[9]={ 0.39894228,0.01328592,0.00225319,
        -0.00157565,0.00916281,-0.02057706,
        0.02635537,-0.01647633,0.00392377};
    float d[9]={ 0.39894228,-0.03988024,-0.00362018,
        0.00163801,-0.01031555,0.02282967,
        -0.02895312,0.01787654,-0.00420059};
    if(n<0) n=-n;
    t=fabs(x);
    if(n!=1)
    {
        if(t<3.75)
        {
            y=(x/3.75)*(x/3.75);
            p=a[6];
            for(i=5; i>=0; i--)	p=p*y+a[i];
        }
        else
        {
            y=3.75/t;
            p=c[8];
            for(i=7; i>=0; i--)	p=p*y+c[i];
            p=p*exp(t)/sqrt(t);
        }
    }
    if(n==0) return(p);
    q=p;
    if(t<3.75)
    {
        y=(x/3.75)*(x/3.75);
        p=b[6];
        for(i=5; i>=0; i--) p=p*y+b[i];
        p=p*t;
    }
    else
    {
        y=3.75/t;
        p=d[8];
        for(i=7; i>=0; i--) p=p*y+d[i];
        p=p*exp(t)/sqrt(t);
    }
    if(x<0.0) p=-p;
    if(n==1) return(p);
    if(doubleEqual(x,0.0)) return(0.0);
    y=2.0/t;
    t = b0 = 0.0;
    b1=1.0;
    m=n+(int)sqrt(40.0*n);
    m=2*m;
    for(i=m; i>0; i--)
    {
        p=b0+i*y*b1;
        b0=b1;
        b1=p;
        if(fabs(b1)>1.0e+10)
        {
            t=t*1.0e-10;
            b0=b0*1.0e-10;
            b1=b1*1.0e-10;
        }
        if(i==n) t=b0;
    }
    p=t*q/b1;
    if((x<0.0)&&(n%2==1)) p=-p;
    return(p);
}
float CalcMethod::TransformativeIntegerBessel2ndFunction(int n, float x)
{
    int i;
    float y,p,b0,b1;

    float a[7]={ -0.57721566,0.4227842,0.23069756,
        0.0348859,0.00262698,0.0001075,0.0000074};
    float b[7]={ 1.0,0.15443144,-0.67278579,
        -0.18156897,-0.01919402,-0.00110404,-0.00004686};
    float c[7]={ 1.25331414,-0.07832358,0.02189568,
        -0.01062446,0.00587872,-0.0025154,0.00053208};
    float d[7]={ 1.25331414,0.23498619,-0.0365562,
        0.01504268,-0.00780353,0.00325614,-0.00068245};
    if(n<0) n=-n;
    if(x<0.0) x=-x;
    if(doubleEqual(x,0.0)) return(1.0e+70);
    if(n!=1)
    {
        if(x < 2.0 || doubleEqual(x, 2.0))
        {
            y=x*x/4.0; p=a[6];
            for(i=5; i>=0; i--) p=p*y+a[i];
            p=p-TransformativeIntegerBessel1stFunction(0,x)*log(x/2.0);
        }
        else
        {
            y=2.0/x; p=c[6];
            for(i=5; i>=0; i--) p=p*y+c[i];
            p=p*exp(-x)/sqrt(x);
        }
    }
    if(n==0) return(p);
    b0=p;
    if(x < 2.0 || doubleEqual(x, 2.0))
    {
        y=x*x/4.0;
        p=b[6];
        for(i=5; i>=0; i--) p=p*y+b[i];
        p=p/x+TransformativeIntegerBessel1stFunction(1,x)*log(x/2.0);
    }
    else
    {
        y=2.0/x;
        p=d[6];
        for(i=5; i>=0; i--) p=p*y+d[i];
        p=p*exp(-x)/sqrt(x);
    }
    if(n==1) return(p);
    b1=p;
    y=2.0/x;
    for(i=1; i<n; i++)
    {
        p=b0+i*y*b1;
        b0=b1;
        b1=p;
    }
    return(p);
}
//...............
//.............................................................................................
//...............
void CalcMethod::waveseparate( complex<float> h[3],float s[3],float SPACE_INTERVAL,float w,complex<float> W[8])
{
    int i,j;
    int NumRows,NumColumns;
    float matrixA[4*8*3],matrixB[2*8];
    float mA[6];
    for(i=0;i<8;i++)
    {
        for(j=0;j<3;j++)
        {
            matrixA[i*6+j]=cos(i*w*s[j]*SPACE_INTERVAL);
        }
        for(j=3;j<6;j++)
        {
            matrixA[i*6+j]=-sin(i*w*s[j-3]*SPACE_INTERVAL);
        }

    }
    for(i=8;i<2*8;i++)
    {
        for(j=0;j<3;j++)
        {
            matrixA[i*6+j]=sin(i*w*s[j]*SPACE_INTERVAL);
        }
        for(j=3;j<6;j++)
        {
            matrixA[i*6+j]=cos(i*w*s[j-3]*SPACE_INTERVAL);
        }

    }
    for(i=0;i<8;i++)
    {
        matrixB[i]=real(W[i]);
    }
    for(i=8;i<16;i++)
    {
        matrixB[i]=imag(W[i-8]);
    }
    NumRows=2*8;NumColumns=2*3;
    CMatrix1 mtxA(NumRows,NumColumns,matrixA);
    NumColumns=1;
    CMatrix1 mtxB(NumColumns,NumColumns,matrixB);
    CLEquations leqs(mtxA,mtxB);
    CMatrix1 mtxResult,mtxQ,mtxR;

    leqs.GetRootsetMqr(mtxResult,mtxQ,mtxR);

    for(i=0;i<3;i++)
    {
        h[i]=complex<float>(mtxResult.GetData()[i],mtxResult.GetData()[i+3]);
    }
}

void CalcMethod::waveseparate2( complex<float> h[2],float s[2],float SPACE_INTERVAL,float w,complex<float> W[8])
{
    int i,j;
    int NumRows,NumColumns;
    float matrixA[4*8*2],matrixB[2*8];
    //         float mA[4];
    for(i=0;i<8;i++)
    {
        for(j=0;j<2;j++)
        {
            matrixA[i*4+j]=cos(i*w*s[j]*SPACE_INTERVAL);
        }
        for(j=2;j<4;j++)
        {
            matrixA[i*4+j]=-sin(i*w*s[j-2]*SPACE_INTERVAL);
        }

    }
    for(i=8;i<2*8;i++)
    {
        for(j=0;j<2;j++)
        {
            matrixA[i*4+j]=sin(i*w*s[j]*SPACE_INTERVAL);
        }
        for(j=2;j<4;j++)
        {
            matrixA[i*4+j]=cos(i*w*s[j-2]*SPACE_INTERVAL);
        }

    }
    for(i=0;i<8;i++)
    {
        matrixB[i]=real(W[i]);
    }
    for(i=8;i<16;i++)
    {
        matrixB[i]=imag(W[i-8]);
    }
    NumRows=2*8;NumColumns=2*2;
    CMatrix1 mtxA(NumRows,NumColumns,matrixA);
    /*float vvv[16*4],voo[16][4],voo1[16][4];
    for(int i=0;i<NumRows*NumColumns;i++)
    {
            vvv[i]=i;
    }
    memcpy(&voo[0],&vvv[0],16*4*4);
    CMatrix1 mtxMM(NumRows,NumColumns,vvv);
    for(int i=0;i<NumRows;i++)
    {
        for(int j=0;j<NumColumns;j++)
        {
            voo1[i][j]=mtxMM.GetElement(i,j);
        }
    }*/
    NumColumns=1;
    CMatrix1 mtxB(NumRows,NumColumns,matrixB);
    CLEquations leqs(mtxA,mtxB);
    CMatrix1 mtxResult,mtxQ,mtxR;

    leqs.GetRootsetMqr(mtxResult,mtxQ,mtxR);

    for(i=0;i<2;i++)
    {
        h[i]=complex<float>(mtxResult.GetData()[i],mtxResult.GetData()[i+2]);
    }
}
//
///三波分离
//
void CalcMethod::WAVEseparate( float wf[8][1024],float hwf1[1024],float hwf2[1024],float hwf3[1024],float dtst,float dfinterval,float SPACE_INTERVAL,int NSAMPLE,float fmin,float fmax)
{
    complex<float> fwf[8][1024],fhf[3],wff[8],fhwf[3][1024],fhwf0[3][1024],temp;
    float xreal[1024],ximag[1024],svp[3];
    int i,j,NFFT=1024;
    int ifmin,ifmax,i0,i1;
    double err,err0,ss,w;
    err0=1.0e18;
    i0=0;
    i1=1;
    ifmin=fmin/dfinterval;
    ifmax=fmax/dfinterval;
    //
    for(i=0;i<8;i++)
    {
        for(j=0;j<NFFT;j++)
        {
            xreal[j]=0.0;
            ximag[j]=0.0;
        }
        for(j=0;j<NSAMPLE;j++)
        {
            xreal[j]=wf[i][j];
        }
        //fft842(int INN, int N, float *X, float *Y)
        fft842(i0,NFFT,xreal,ximag);
        for(j=0;j<NFFT;j++)
        {
            fwf[i][j]=complex<float>(xreal[j],ximag[j]);
        }
    }
    //
    for( ss=dtst-10;ss<dtst+10;ss=ss+2.)
    {
        svp[0]=0.0;svp[1]=2*ss/1000000;svp[2]=-2*ss/100000;
        err=0.0;
        for(i=ifmin;i<ifmax;i++)
        {
            w=2*3.1415926*(i-1)*dfinterval;
            for(j=0;j<8;j++)wff[j]=fwf[j][i];//complex<float>(sqrt(real(fwf[j][i])*real(fwf[j][i])+real(fwf[j][i])*real(fwf[j][i])),0.0);
            //waveseparate(  h[3],float s[3],float SPACE_INTERVAL,float w,complex<float> W[8])
            waveseparate3( fhf,svp, SPACE_INTERVAL, w,wff);
            for(j=0;j<3;j++)fhwf0[j][i]=fhf[j];
            temp=fhf[0]-fhf[1]-fhf[2];
            float temp1=sqrt(real(temp)*real(temp)+imag(temp)*imag(temp));
            float temp2=sqrt(real(wff[0])*real(wff[0])+imag(wff[0])*imag(wff[0]));
            err=err+(temp1-temp2)*(temp1-temp2);
        }
        if(err<err0)
        {
            err0=err;
            for(i=ifmin;i<ifmax;i++)
            {
                for(j=0;j<3;j++)fhwf[j][i]=fhwf0[j][i];
            }
        }

    }

    for(i=0;i<3;i++)
    {
        for(j=0;j<NFFT;j++)
        {
            xreal[j]=0.0;//real(fwf[0][j]);
            ximag[j]=0.0;//imag(fwf[0][j]);
        }
        for(j=ifmin;j<ifmax;j++)
        {
            xreal[j]=real(fhwf[i][j]);
            ximag[j]=imag(fhwf[i][j]);
            //            xreal[NFFT-j+1]=real(fhwf[i][j]);
            //			ximag[NFFT-j+1]=imag(fhwf[i][j]);
        }
        //fft842(int INN, int N, float *X, float *Y)
        fft842(i1,NFFT,xreal,ximag);
        for(j=0;j<NSAMPLE;j++)
        {
            if(i==0)	hwf1[j]=xreal[j];
            if(i==1)	hwf2[j]=xreal[j];
            if(i==2)	hwf3[j]=xreal[j];
        }
    }

}
//
///二波分离
//
void CalcMethod::WAVEseparate2( float wf[8][1024],float hwf1[1024],float hwf2[1024],float dtst,float dfinterval,float SPACE_INTERVAL,int NSAMPLE,int iflag,float fmin,float fmax)
{
    complex<float> fwf[8][1024],fhf[2],wff[8],fhwf[2][1024],fhwf0[2][1024],temp;
    float xreal[1024],ximag[1024],svp[2];
    int i,j,NFFT=1024;
    int ifmin,ifmax,i0,i1;
    double err,err0,ss,w;
    err0=1.0e20;
    i0=0;
    i1=1;
    ifmin=fmin/dfinterval;
    ifmax=fmax/dfinterval;
    //  ifmin=0;
    //	ifmax=NFFT/2;
    //
    for(i=0;i<8;i++)
    {
        for(j=0;j<NFFT;j++)
        {
            xreal[j]=0.0;
            ximag[j]=0.0;
        }
        for(j=0;j<NSAMPLE;j++)
        {
            xreal[j]=wf[i][j];
        }
        //fft842(int INN, int N, float *X, float *Y)
        fft842(i0,NFFT,xreal,ximag);
        for(j=0;j<NFFT;j++)
        {
            fwf[i][j]=complex<float>(xreal[j],ximag[j]);
        }
    }
    //
    for( ss=dtst-15;ss<dtst+15;ss=ss+2.5)
    {
        if(iflag==0)
        {
            svp[0]=-ss/1000000;svp[1]=ss/1000000;
        }
        if(iflag==1)
        {
            svp[0]=0;svp[1]=2*ss/1000000;
        }
        if(iflag==2)
        {
            svp[0]=2*ss/1000000;svp[1]=-2*ss/1000000;
        }
        err=0.0;

        for(i=ifmin;i<ifmax;i++)
        {
            w=2*3.1415926*(i)*dfinterval;
            for(j=0;j<8;j++)wff[j]=fwf[j][i];
            //complex<float>(sqrt(real(fwf[j][i])*real(fwf[j][i])+real(fwf[j][i])*real(fwf[j][i])),0.0);
            //waveseparate(  h[3],float s[3],float SPACE_INTERVAL,float w,complex<float> W[8])
            //wavesp2(  fhf,svp, SPACE_INTERVAL, w,wff);
            waveseparate2(  fhf,svp, SPACE_INTERVAL, w,wff);
            for(j=0;j<2;j++)fhwf0[j][i]=fhf[j];
            float temp0=sqrt(real(wff[0]-fhf[0]-fhf[1])*real(wff[0]-fhf[0]-fhf[1])+imag(wff[0]-fhf[0]-fhf[1])*imag(wff[0]-fhf[0]-fhf[1]));
            //          float temp0=sqrt(real(fhf[0])*real(fhf[0])+imag(fhf[0])*imag(fhf[0]));
            //		   float temp1=sqrt(real(fhf[1])*real(fhf[1])+imag(fhf[1])*imag(fhf[1]));
            //          float temp2=sqrt(real(wff[0])*real(wff[0])+imag(wff[0])*imag(wff[0]));
            err=err+temp0;//fabs(temp2-temp1-temp0);//sqrt(real(temp)*real(temp)+imag(temp)*imag(temp));
        }
        if(err<err0)
        {
            err0=err;
            for(i=ifmin;i<ifmax;i++)
            {
                for(j=0;j<2;j++)fhwf[j][i]=fhwf0[j][i];
            }
        }


    }

    for(i=0;i<2;i++)
    {
        for(j=0;j<NFFT;j++)
        {
            xreal[j]=0.0;//real(fwf[0][j]);
            ximag[j]=0.0;//imag(fwf[0][j]);
        }
        for(j=ifmin;j<ifmax;j++)
        {
            xreal[j]=real(fhwf[i][j]);
            ximag[j]=imag(fhwf[i][j]);

            xreal[NFFT-j-1]=real(fhwf[i][j]);
            ximag[NFFT-j-1]=-imag(fhwf[i][j]);

        }
        //fft842(int INN, int N, float *X, float *Y)
        fft842(i1,NFFT,xreal,ximag);
        for(j=0;j<NSAMPLE;j++)
        {
            if(i==0)	hwf1[j]=xreal[j];
            if(i==1)	hwf2[j]=xreal[j];
            //		   if(i==2)	hwf3[j]=xreal[j];
        }
    }

}

void CalcMethod::waveseparate3(complex<float> h[3],float s[3],float SPACE_INTERVAL,float w,complex<float> W[8])
{
    complex<float> Q[8][8], Z[8][3];
    int i,j;
    for(i=0;i<8;i++)
    {
        for(j=0;j<3;j++)
        {
            Z[i][j]=complex<float>(cos(i*w*SPACE_INTERVAL*s[j]),sin(i*w*SPACE_INTERVAL*s[j]));
        }
    }
    splitCQR(Z, Q);
    setrootQR(Z,h,Q,W);

}

void  CalcMethod::setrootQR(complex<float> R[8][3],complex<float> h[],complex<float> Q[8][8],complex<float> B[])
{
    complex<float> R1[8][3],QT[8][8],C[8];
    int i,j;
    for(i=0;i<3;i++)
    {
        for(j=0;j<3;j++)
        {
            R1[i][j]=R[i][j];
        }
    }
    for(i=0;i<8;i++)
    {
        for(j=0;j<8;j++)
        {
            QT[i][j]=conj(Q[j][i]);
        }
    }
    for(i=0;i<8;i++)
    {
        C[i]=0;
        for(j=0;j<8;j++)
        {
            C[i]=C[i]+QT[i][j]*B[j];
        }
    }
    h[2]=C[2]/R1[2][2];
    h[1]=(C[i]-h[2]*R1[1][2])/R1[1][1];
    h[0]=(C[i]-h[2]*R1[0][2]-h[1]*R1[0][1])/R1[0][0];



}

void  CalcMethod::splitCQR(complex<float> Z[8][3],complex<float> Q[8][8])
{
    int i,j,p,q,m;
    complex<float> z[3][3],h1[3][3],h[8][8],h2[8][8];
    for(i=0;i<3;i++)
    {
        for(j=0;j<3;j++)
        {
            z[i][j]=Z[i][j];
        }
    }
    splitC3QR(z,h1);
    for(i=0;i<8;i++)
    {
        for(j=0;j<8;j++)
        {
            if(i<3&&j<3)
            {
                Q[i][j]=h1[i][j];
            }
            else
            {
                Q[i][j]=Z[i][j];
            }
        }
    }
    for(i=0;i<3;i++)
    {
        for(j=0;j<3;j++)
        {
            Z[i][j]=z[i][j];
        }
    }
    for(i=3;i<8;i++)
    {
        recursion(Z,i,h);
        for(p=0;p<8;p++)
        {
            for(q=0;q<8;q++)
            {

                for(m=0;m<8;m++)
                {
                    h2[p][q]=complex<float>(0,0);
                    h2[p][q]=h2[p][q]+h[p][m]*Q[m][q];
                }
            }
        }
        for(p=0;p<8;p++)
        {
            for(q=0;q<8;q++)
            {
                Q[p][q]=h2[p][q];
            }
        }
    }


}
void  CalcMethod::recursion(complex<float> Z[8][3],int n,complex<float> Q[8][8])
{
    float e[3],o[3];
    complex<float> I[8][8],q[3],q1[3],u[3],r[3],h[8][8],h1[8][8],u1[8];
    int i,j,k,l,m;
    for(i=0;i<8;i++)
    {
        for(j=0;j<8;j++)
        {
            if(i==j)
            {
                I[i][j]=complex<float>(1,0);
                h1[i][j]=complex<float>(1,0);
            }
            else
            {
                I[i][j]=complex<float>(0,0);
                h1[i][j]=complex<float>(0,0);
            }
        }
    }
    for(k=0;k<3;k++)
    {
        e[k]=sqrt(pow(abs(Z[k][k]),2)+pow(abs(Z[n][k]),2));
        o[k]=e[k]*(e[k]+abs(Z[k][k]));
        q[k]=Z[k][k]+complex<float>(e[k]*cos(arg(Z[k][k])),e[k]*sin(arg(Z[k][k])));
        q1[k]=q[k]/complex<float>(o[k],0);
        u[k]=Z[n][k]/complex<float>(o[k],0);
        for(l=k;l<4;l++)
        {
            r[l]=conj(q[k])*Z[k][l]+conj(Z[n][k])*Z[n][l];
            Z[k][l]=Z[k][l]-q1[k]*r[l];
            Z[n][l]=Z[n][l]-u[k]*r[l];
        }
        for(i=0;i<8;i++)
        {
            if(i==k)
            {
                u1[i]=q[k];
            }
            else if(i==n)
            {
                u1[i]=u[k];
            }
            else
            {
                u1[i]=complex<float>(0,0);
            }

        }
        for(i=0;i<8;i++)
        {
            for(j=0;j<8;j++)
            {
                h[i][j]=I[i][j]-u1[i]*conj(u1[j]);
            }
        }
        for(i=0;i<8;i++)
        {
            for(j=0;j<8;j++)
            {
                for(m=0;m<8;m++)
                {
                    Q[i][j]=complex<float>(0,0);
                    Q[i][j]=Q[i][j]+h[i][m]*h1[m][j];
                }
            }
        }
        for(i=0;i<8;i++)
        {
            for(j=0;j<8;j++)
            {
                h1[i][j]=h[i][j];
            }
        }


    }
}
void  CalcMethod::splitC3QR(complex<float> z[3][3],complex<float> h[3][3])
{
    float p;
    complex<float> u1[3],u2[2],u3[3];
    complex<float> H[3][3],H1[2][2],H2[3][3];
    complex<float> I[3][3],I2[2][2];
    complex<float> z1[3][3],z2[2][2],z3[2][2];
    complex<float> s[3][3];
    p=sqrt(pow(abs(z[0][0]),2)+pow(abs(z[1][0]),2)+pow(abs(z[2][0]),2));
    u1[0]=z[0][0]+complex<float>(p*cos(arg(z[0][0])),p*sin(arg(z[0][0])));
    u1[1]=z[1][0];
    u1[2]=z[2][0];
    int i,j,m;
    for(i=0;i<3;i++)
    {
        for(j=0;j<3;j++)
        {
            if(i==j)
            {
                I[i][j]=complex<float>(1,0);
            }
            else
            {
                I[i][j]=complex<float>(0,0);
            }
        }
    }
    for(i=0;i<3;i++)
    {
        for(j=0;j<3;j++)
        {
            H[i][j]=I[i][j]-u1[i]*conj(u1[j])/complex<float>(p*(p+abs(z[0][0])),0);
        }
    }
    for(i=0;i<3;i++)
    {
        for(j=0;j<3;j++)
        {
            z1[i][j]=H[i][0]*z[0][j]+H[i][1]*z[1][j]+H[i][2]*z[2][j];
        }
    }
    for(i=0;i<3;i++)
    {
        for(j=0;j<3;j++)
        {
            z[i][j]=z1[i][j];
        }
    }
    z2[0][0]=z[1][1];
    z2[0][1]=z[1][2];
    z2[1][0]=z[2][1];
    z2[0][1]=z[2][2];

    p=sqrt(pow(abs(z1[1][1]),2)+pow(abs(z1[2][1]),2));
    u2[0]=z1[1][1]+complex<float>(p*cos(arg(z1[1][1])),p*sin(arg(z1[2][1])));
    u2[1]=z1[2][1];
    I2[0][0]=complex<float>(1,0);
    I2[0][1]=complex<float>(0,0);
    I2[1][0]=complex<float>(0,0);
    I2[1][1]=complex<float>(1,0);
    for(i=0;i<2;i++)
    {
        for(j=0;j<2;j++)
        {
            H1[i][j]=I2[i][j]-u2[i]*conj(u2[j])/complex<float>(p*(p+abs(z1[1][1])),0.0);
        }
    }
    for(i=0;i<2;i++)
    {
        for(j=0;j<2;j++)
        {
            z3[i][j]=H1[i][0]*z2[0][j]+H1[i][1]*z2[1][j];
        }
    }
    z[1][1]=z3[0][0];
    z[1][2]=z3[0][1];
    z[2][1]=z3[1][0];
    z[2][2]=z3[1][1];
    h[0][0]=complex<float>(1,0);
    h[0][1]=complex<float>(0,0);
    h[0][2]=complex<float>(0,0);
    h[1][0]=complex<float>(0,0);
    h[2][0]=complex<float>(0,0);
    h[1][1]=H1[0][0];
    h[1][2]=H1[0][1];
    h[2][1]=H1[1][0];
    h[2][2]=H1[1][1];
    for(i=0;i<3;i++)
    {
        for(j=0;j<3;j++)
        {
            H2[i][j]=complex<float>(0,0);
            for(m=0;m<3;m++)
            {
                H2[i][j]=H2[i][j]+h[i][m]*H[m][j];
            }
        }
    }
    for(i=0;i<3;i++)
    {
        for(j=0;j<3;j++)
        {
            h[i][j]=H2[i][j];
        }
    }
}

void CalcMethod::wavesp2(complex<float> h[2],float s[2],float d,float w,complex<float> W[8])
{
    complex<float> Z[8][2];
    complex<float> Q[8][8];
    int i,j;
    for(i=0;i<8;i++)
    {
        for(j=0;j<2;j++)
        {
            Z[i][j]=complex<float>(cos(i*w*s[j]*d),sin(i*w*s[j]*d));
        }
    }
    splitQR2(Z,Q);
    setrootQR2(Z,Q,h,W);

}
void CalcMethod::setrootQR2(complex<float> Z[8][2],complex<float> Q[8][8], complex<float> h[2],complex<float> W[8])
{
    int i,j;
    complex<float> C[8];
    for(i=0;i<8;i++)
    {
        C[i]=0;
        for(j=0;j<8;j++)
        {
            C[i]=C[i]+Q[i][j]*W[j];
        }
    }
    h[1]=C[1]/Z[1][1];
    h[0]=(C[0]-Z[0][1]*h[1])/Z[0][0];

}
void CalcMethod::splitQR2(complex<float> Z[8][2],complex<float> Q[8][8])
{
    int i,j,m;
    float p=0.0;
    complex<float> q1,q2;
    complex<float> u1[8],u2[8];
    complex<float> H[8][8],H1[8][8],I[8][8],Z1[8][2],Z2[8][2];
    for(i=0;i<8;i++)
    {
        for(j=0;j<8;j++)
        {
            if(i==j)
            {
                I[i][j]=complex<float>(1,0);
            }
            else
            {
                I[i][j]=complex<float>(0,0);
            }

        }
    }
    for(i=0;i<8;i++)
    {
        p=p+pow(abs(Z[i][0]),2);
    }
    p=sqrt(p);

    q1=complex<float>(p*cos(arg(Z[0][0])),p*sin(arg(Z[0][0])));
    u1[0]=Z[0][0]+complex<float>(p*cos(arg(Z[0][0])),p*sin(arg(Z[0][0])));

    for(i=1;i<8;i++)
    {
        u1[i]=Z[i][0];
    }
    for(i=0;i<8;i++)
    {
        for(j=0;j<8;j++)
        {
            H[i][j]=I[i][j]-u1[i]*conj(u1[j])/(q1*(q1+Z[0][0]));
        }
    }
    for(i=0;i<8;i++)
    {
        for(j=0;j<2;j++)
        {
            Z1[i][j]=complex<float>(0,0);
            for(m=0;m<8;m++)
            {
                Z1[i][j]=Z1[i][j]+H[i][m]*Z[m][j];
            }
        }
    }
    p=0;
    for(i=1;i<8;i++)
    {
        p=p+pow(abs(Z1[i][1]),2);
    }
    q2=complex<float>(p*cos(arg(Z[1][1])),p*sin(arg(Z[1][1])));
    u2[0]=complex<float>(0,0);
    u2[1]=Z1[1][1]+q2;
    for(i=2;i<8;i++)
    {
        u2[i]=Z[i][1];
    }
    for(i=0;i<8;i++)
    {
        for(j=0;j<8;j++)
        {
            H1[i][j]=I[i][j]-u2[i]*conj(u2[j])/(q2*(q2+Z[1][1]));
        }
    }
    for(i=0;i<8;i++)
    {
        for(j=0;j<2;j++)
        {
            Z2[i][j]=complex<float>(0,0);
            for(m=0;m<8;m++)
            {
                Z2[i][j]=Z2[i][j]+H1[i][m]*Z1[m][j];
            }
        }
    }
    for(i=0;i<8;i++)
    {
        for(j=0;j<2;j++)
        {
            Z[i][j]=Z2[i][j];
        }
    }
    for(i=0;i<8;i++)
    {
        for(j=0;j<8;j++)
        {
            Q[i][j]=complex<float>(0,0);
            for(m=0;m<8;m++)
            {
                Q[i][j]=Q[i][j]+H1[i][m]*H[m][j];
            }
        }
    }

}

float CalcMethod::Func1(float x,float y ,float Vpv,float Vsv,float Vsh,float vf,float denf,
                        float den,float R,float w,int n,float MTOOL,float R_TOOL)
{
    float c11,c13,c33,c44,c66;
    float etool,temp;
    complex<float> k;
    complex<float> U,V,W;
    complex<float> Qp,Qsh,Qsv;
    complex<float> Q[4][4];
    complex<float> fun,T,T1,T2;
    complex<float> _c11,_c13,_c66,_c44,_c33;
    complex<float> _R,_n,_denf,_w,_den,_vf;
    complex<float> _a,_b,f,kf;
    complex<float> a,ai,_ar,ttt;
    int t,i,j,i1,i0;
    int u[24][4]={{0,1,2,3},{0,1,3,2},{0,2,1,3},{0,2,3,1},{0,3,1,2},{0,3,2,1},
    {1,0,2,3},{1,0,3,2},{1,2,0,3},{1,2,3,0},{1,3,0,2},{1,3,2,0},
    {2,0,3,1},{2,0,1,3},{2,1,3,0},{2,1,0,3},{2,3,0,1},{2,3,1,0},
    {3,0,1,2},{3,0,2,1},{3,1,2,0},{3,1,0,2},{3,2,0,1},{3,2,1,0}};

    k=complex<float>(x,y);
    i1=1;i0=0;
    c44=den*Vsv*Vsv;
    c66=den*Vsh*Vsh;
    c33=den*Vpv*Vpv;
    c11=(c66/c44)*c33;
    c13=c33-2*c44;


    a=complex<float>(R_TOOL,0.0);
    ai=complex<float>(0.0,1.0);
    ai=complex<float>(-1.0,0.0);
    ttt=complex<float>(1.0e-10,0.0);
    temp=MTOOL/R_TOOL;

    _c11=complex<float>(c11,0);
    _c13=complex<float>(c13,0);
    _c66=complex<float>(c66,0);
    _c44=complex<float>(c44,0);
    _c33=complex<float>(c33,0);
    _R=complex<float>(R,0);
    _n=complex<float>(n,0);
    _denf=complex<float>(denf,0);
    _w=complex<float>(w,0);
    _den=complex<float>(den,0);
    _vf=complex<float>(vf,0);

    U=_c11*_c44;
    V=_den*(_c11+_c44)-(_c11*_c33-_c13*_c13-complex<float>(2.0,0.0)*_c13*_c44)*(k*k/(_w*_w));
    W=_c33*_c44*(_den/_c44-k*k/(_w*_w))*(_den/_c33-k*k/(_w*_w));

    T=(_c44*k*k-_den*_w*_w)/_c66;
    Qsh=sqrt(T);
    T1=V*V-complex<float>(4.0,0.0)*U*W;
    T2=(complex<float>(-1.0,0.0)*V+sqrt(T1))/(complex<float>(2.0,0.0)*U);
    Qp=_w*sqrt(T2);
    T2=(complex<float>(-1.0,0.0)*V-sqrt(T1))/(complex<float>(2.0,0.0)*U);
    Qsv=_w*sqrt(T2);

    _a=(complex<float>(-1.0,0.0)/(ai*k))*((_c13+complex<float>(2.0,0.0)*_c44)*k*k-
        _c11*Qp*Qp-_den*_w*_w)/(_c44*k*k-(_c11-_c13-_c44)*Qp*Qp-_den*_w*_w);
    _b=_ar*ai*k*(_c44*k*k-(_c11-_c13-_c44)*Qsv*Qsv-_den*_w*_w)/
        ((_c13+complex<float>(2.0,0.0)*_c44)*k*k-_c11*Qsv*Qsv-_den*_w*_w);
    kf=_w/_vf;
    T=k*k-kf*kf;
    f=sqrt(T);

    Q[0][0]=_ar*(_n/_R)*Ibessel(n,f*_R)-f*Ibessel(n+1,f*_R);
    Q[0][1]=_ar*(complex<float>(1.0,0.0)+ai*k*_a)*Qp*((_ar*_n/(Qp*_R))*Kbessel(n,Qp*_R)+Kbessel(n+1,Qp*_R));
    Q[0][2]=(_n/_R)*Kbessel(n,Qsh*_R);
    Q[0][3]=_ar*(ai*k+_b)*Qsv*((_ar*_n/(Qp*_R))*Kbessel(n,Qp*_R)+Kbessel(n+1,Qp*_R));
    Q[1][0]=_denf*_w*_w*Ibessel(n,f*_R);
    Q[1][1]=(_c11*Qp*Qp-_c13*k*k+(_c11-_c13)*ai*k*_a*Qp*Qp)*Kbessel(n,Qp*_R)+complex<float>(2.0,0.0)*_c66*(Qp/_R)*(complex<float>(1.0,0.0)+ai*k*_a)*((complex<float>(n*(n-1),0.0)/(Qp*_R))*Kbessel(n,Qp*_R)+Kbessel(n+1,Qp*_R));
    Q[1][2]=complex<float>(-2,0)*_c66*_n*(Qsh/_R)*((complex<float>(1-n,0.0)/(Qsh*_R))*Kbessel(n,Qsh*_R)+Kbessel(n+1,Qsh*_R));
    Q[1][3]=((_c11*Qsv*Qsv-_c13*k*k)*_b+(_c11-_c13)*ai*k*Qsv*Qsv)*Kbessel(n,Qsv*_R)+complex<float>(2,0)*_c66*(Qsv/_R)*(ai*k+_b)*((complex<float>(n*(n-1),0.0)/(Qsv*_R))*Kbessel(n,Qsv*_R)+Kbessel(n+1,Qsv*_R));
    Q[2][0]=complex<float>(0,0);
    Q[2][1]=complex<float>(2,0)*_c66*_n*(Qp/_R)*(complex<float>(1.0,0.0)+ai*k*_a)*((complex<float>(1-n,0.0)/(Qp*_R))*Kbessel(n,Qp*_R)+Kbessel(n,Qp*_R));
    Q[2][2]=_ar*_c66*Qsh*Qsh*((complex<float>(1.0,0.0)+complex<float>(2.0*n*(n-1),0.0)/(Qsh*_R*Qsh*_R))*Kbessel(n,Qsh*_R)+(complex<float>(2.0,0.0)/(Qsh*_R))*Kbessel(n+1,Qsh*_R));
    Q[2][3]=complex<float>(2,0)*_c66*_n*(Qsv/_R)*(ai*k+_b)*((complex<float>(1-n,0.0)/(Qsv*_R))*Kbessel(n,Qsv*_R)+Kbessel(n+1,Qsv*_R));
    Q[3][0]=complex<float>(0,0);
    Q[3][1]=_ar*_c44*Qp*(complex<float>(2,0)*ai*k-_a*(k*k+Qp*Qp))*(_ar*(_n/(Qp*_R))*Kbessel(n,Qp*_R)+Kbessel(n+1,Qp*_R));
    Q[3][2]=complex<float>(1,0)*ai*k*_n*(_c44/_R)*Kbessel(n,Qsh*_R);
    Q[3][3]=_c44*Qsv*((k*k+Qsv*Qsv)-complex<float>(2,0)*ai*k*_b)*((_ar*_n/(Qsv*_R))*Kbessel(n,Qsv*_R)+Kbessel(n+1,Qsv*_R));//Kbessel(1,Qsv*_R);

    etool=(temp*real(f*Ibessel(i1,f*a))+denf*w*w*real(Ibessel(i0,f*a)))/(temp*real(f*Kbessel(i1,f*a))-denf*w*w*real(Kbessel(i0,f*a)));
    if(MTOOL>0.01)
    {
        if(n==0)
        {
            Q[0][0]=_ar*(_n/_R)*Ibessel(n,f*_R)-f*(Ibessel(n+1,f*_R)-complex <float>(etool,0)*Kbessel(n+1,f*_R));
            Q[1][0]=_denf*_w*_w*(Ibessel(n,f*_R)+complex <float>(etool,0)*Kbessel(n,f*_R));
        }
        if(n==1)
        {
            Q[0][0]=_ar*(_n/_R)*(Ibessel(n,f*_R)-complex <float>(etool,0)*Kbessel(n,f*_R))-f*(Ibessel(n+1,f*_R));
            Q[1][0]=_denf*_w*_w*(Ibessel(n,f*_R)-complex <float>(etool,0)*Kbessel(n,f*_R));
        }

    }
    for(j=1;j<=3;j++)
    {
        for(i=0;i<=3;i++)Q[j][i]=Q[j][i]*ttt;
    }


    for(j=0;j<24;j++)
    {
        t=reversion(u[j]);
        fun=fun+complex<float>(pow(-1.0,t),0)*Q[0][u[j][0]]*Q[1][u[j][1]]*Q[2][u[j][2]]*Q[3][u[j][3]];

    }
    return sqrt(real(fun)*real(fun)+imag(fun)*imag(fun));

}


complex<float> CalcMethod::GetRootMonteCarlo1(float Vpv,float Vsv,float Vsh,float vf,float denf,float den,float R,float w,int n,float MTOOL,float R_TOOL)
{
    int k;
    float xx,yy,a,r,z,x1,y1,z1,eps=1.0e-10,k_real,k_imag;
    int nControlB=10;
    float xStart;
    //	float PI=3.1415926;
    complex<float> kw0,kw;

    //	w=2*PI*fff;

    if(Vsv>vf)
    {
        kw0=complex <float>(w/vf+.10,0.0);
    }
    else
    {
        kw0=complex <float>(w/vf-.10,0.0);
    }
    xStart=0.10;//w/vf;
    if(n==1)
    {
        kw0=complex <float>(w/vf-.50,0.0);
        xStart=0.50;//w/vf;
    }
    a=xStart;
    k=1;
    r=1.0;
    xx=real(kw0);
    yy=imag(kw0);
    z=Func1(xx,yy,Vpv,Vsv,Vsh,vf,denf,den,R, w,n,MTOOL,R_TOOL);

    // ???è?????ó?vp
    while (a>=eps)
    {
        x1=-a+2.0*a*rnd(&r);
        x1=xx+x1;
        y1=-a+2.0*a*rnd(&r);
        y1=yy+y1;

        z1=Func1(x1,y1,Vpv,Vsv,Vsh,vf,denf,den,R,w,n,MTOOL,R_TOOL);

        k=k+1;
        if (z1>=z)
        {
            if (k>nControlB)
            {
                k=1;
                a=a/2.0;
            }
        }
        else
        {
            //			k=1;
            xx=x1;
            yy=y1;
            z=z1;
            if (z<eps)
            {
                k_real=xx;
                k_imag=yy;
                kw=complex <float>(xx,yy);
                return kw;
            }
        }
    }

    k_real=xx;
    k_imag=yy;
    kw=complex<float>(xx,yy);
    return kw;
}

complex<float> CalcMethod::xianjie(float vp,float vs,float vf,float denf,float den,float R,float w,int n,float MTOOL,float R_TOOL)
{
    complex<float> xk ;
    complex<float> xk_1;
    complex<float>  x,fk,fk_1;
    float e=0.00001;
    int k=0;
    complex<float> x0, x1;
    if(n==0)
    {
        x0=complex<float>(w/vf/0.99,0.0);
        x1=complex<float>(w/vf/0.9,0.0);
    }
    if(n==1)
    {
        x0=complex<float>(1.1*w/vs,0.0);
        x1=complex<float>(0.9*w/vf,0.0);
    }

    xk=x1,xk_1=x0;
    while(abs(xk-xk_1)>e)
    {
        float temp=abs(xk-xk_1);
        fk=Funcxianjie(real(xk),imag(xk) ,vp,vs,vf,denf,den,R,w,n,MTOOL,R_TOOL);
        fk_1=Funcxianjie(real(xk_1),imag(xk_1),vp,vs,vf,denf,den,R,w,n,MTOOL,R_TOOL);
        x=xk-fk*(xk-xk_1)/(fk-fk_1);
        xk_1=xk;
        xk=x;
        if(k>100)return x;
        k++;
    }
    return x;

}


complex<float> CalcMethod::Funcxianjie(float x,float y ,float vp,float vs,float vf,float denf,float den,float R,float w,int n,float MTOOL,float R_TOOL)
{
    complex<float>  M[4][4],fun,p,s,f,g,k,ttt;
    complex<float> fr,sr,pr,fa,ai,inter[4],Inter;
    complex<float> common;
    float arg0=0.0,SPACE=3.6;
    float kp,kf,ks,etool,temp;
    int j,i;
    //       k=w/v;
    k=complex<float>(x,y);
    kp=w/vp;
    ks=w/vs;
    kf=w/vf;
    ai=complex<float>(0.0,1.0);
    ttt=complex<float>(1.e-10,0.0);
    p=sqrt(k*k-complex<float>(kp*kp,0));
    s=sqrt(k*k-complex<float>(ks*ks,0));
    f=sqrt(k*k-complex<float>(kf*kf,0));

    inter[0]=pow(f*complex<float>(R_TOOL/2,0),n);
    inter[1]=complex<float>(cos(n*(arg0-3.1415926/2)),0);
    Inter=ai*k*complex<float>(SPACE,0);
    inter[2]=complex<float>(pow(2.718282,(double)real(Inter))*cos(imag(Inter)),pow(2.718282,(double)real(Inter))*sin(imag(Inter)));
    inter[3]=inter[0]*inter[1]*inter[2];
    common=complex<float>(1.0,0.0);//inter[3]*inter[3]*inter[3]*inter[3]*complex<float>(1.e-12,0.0);


    fr=f*complex <float>(R,0.0);
    pr=p*complex <float>(R,0.0);
    sr=s*complex <float>(R,0.0);
    fa=f*complex <float>(R_TOOL,0.0);

    temp=MTOOL/R_TOOL;
    etool=(temp*real(f*Ibessel(1,fa))+denf*w*w*real(Ibessel(0,fa)))/(temp*real(f*Kbessel(1,fa))-denf*w*w*real(Kbessel(0,fa)));
    M[0][0]=complex <float>(-1.,0.0)*complex <float>(n/R,0.0)*Ibessel(n,fr)-f*Ibessel(n+1,fr);
    M[0][1]=complex <float>(-1.,0.0)*p*(complex <float>(-n,0)/(pr)*Kbessel(n,pr)+Kbessel(n+1,pr));
    M[0][2]=complex <float>(-1.,0.0)*complex <float>(n/R,0)*Kbessel(n,sr);
    M[0][3]=ai*complex <float>(-ks,0.0)*(complex <float>(-n,0.0)/(sr)*Kbessel(n,sr)+Kbessel(n+1,sr));
    M[1][0]=(complex <float>(denf*w*w,0.0)*Ibessel(n,fr));
    M[1][1]=(complex <float>(den,0)*(complex <float>(2.0*vs*vs,0)*k*k-complex <float>(w*w,0.0))*Kbessel(n,pr)+complex <float>(2.0*den*vs*vs/R,0.0)*p*((complex <float>(n*(n-1),0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[1][2]=(complex <float>(-2*n*den*vs*vs/R,0.0)*s*((complex <float>(1-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    M[1][3]=((ai*k*complex <float>(2*den*vs*vs,0.0)*s*s*Kbessel(n,sr)+ai*complex <float>(2*ks*den*vs*vs/R,0)*((complex <float>(n*(n-1),0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr))));
    M[2][0]=(0.0,0.0);
    M[2][1]=(complex <float>(-2*n*den*vs*vs,0.0)*p*((complex <float>(1-n,0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[2][2]=(complex <float>(-den*vs*vs,0)*s*s*((complex <float>(1.0,0.0)+complex <float>(2*n*(n-1),0.0)/(sr*sr))*Kbessel(n,sr)+(complex <float>(2,0.0)/sr)*Kbessel(n+1,sr)));
    M[2][3]=(ai*complex <float>(2*ks*n*den*vs*vs/R,0)*((complex <float>(1-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    M[3][0]=(0.0,0.0);
    M[3][1]=(ai*k*complex <float>(-2*den*vs*vs,0)*p*((complex <float>(-n,0.0)/pr)*Kbessel(n,pr)+Kbessel(n+1,pr)));
    M[3][2]=(ai*k*complex <float>(-n*den*vs*vs/R,0)*Kbessel(n,sr));
    M[3][3]=((k*k+s*s)*s*complex <float>(den*vs*vs,0)*((complex <float>(-n,0)/sr)*Kbessel(n,sr)+Kbessel(n+1,sr)));
    if(MTOOL>0.01)
    {
        if(n==0)
        {
            M[0][0]=complex <float>(-1.,0.0)*complex <float>(n/R,0.0)*Ibessel(n,fr)-f*(Ibessel(n+1,fr)-etool*Kbessel(n+1,fr));
            M[1][0]=(complex <float>(denf*w*w,0.0)*(Ibessel(n,fr)+etool*Kbessel(n,fr)));
        }
        if(n==1)
        {
            M[0][0]=complex <float>(-1.,0.0)*complex <float>(n/R,0.0)*(Ibessel(n,fr)-etool*Kbessel(n,fr))-f*(Ibessel(n+1,fr));
            M[1][0]=(complex <float>(denf*w*w,0.0)*(Ibessel(n,fr)-etool*Kbessel(n,fr)));
        }

    }
    for(j=0;j<=3;j++)
    {
        for(i=0;i<=3;i++)
        {
            M[j][i]=M[j][i]*ttt;

        }
    }

    int u[24][4]={{0,1,2,3},{0,1,3,2},{0,2,1,3},{0,2,3,1},{0,3,1,2},{0,3,2,1},
    {1,0,2,3},{1,0,3,2},{1,2,0,3},{1,2,3,0},{1,3,0,2},{1,3,2,0},
    {2,0,3,1},{2,0,1,3},{2,1,3,0},{2,1,0,3},{2,3,0,1},{2,3,1,0},
    {3,0,1,2},{3,0,2,1},{3,1,2,0},{3,1,0,2},{3,2,0,1},{3,2,1,0}};
    int t,sd[4];
    for(j=0;j<24;j++)
    {
        sd[0]=u[j][0];
        sd[1]=u[j][1];
        sd[2]=u[j][2];
        sd[3]=u[j][3];

        t=reversion(sd);
        inter[0]=M[0][u[j][0]];
        inter[1]=inter[0]*M[1][u[j][1]];
        inter[2]=inter[1]*M[2][u[j][2]];
        inter[3]=inter[2]*M[3][u[j][3]];
        fun=fun+complex <float>(pow(-1.0,t),0.0)*inter[3]*common;

    }
    return fun;
}
complex<float> CalcMethod::interposek(float x,float y ,float vp,float vs,float vf,float denf,float den,float R,float w,int n,float MTOOL,float R_TOOL)
{
    float hr=0.1;//x/10;//?
    float hi=0.0;//y/10;//?
    complex<float> X[5];
    complex<float> inter[8],Inter;
    int i;
    for(i=0;i<5;i++)
    {
        X[i]=complex<float>(x+(i)*hr,y+(i)*hi);
    }
    complex<float> h=complex<float>(hr,hi);
    complex<float> m,f0,f1,f2,f3,f4;
    f0=Funcxianjie(real(X[0]),imag(X[0]),vp,vs,vf,denf,den,R,w,n,MTOOL,R_TOOL);
    f1=Funcxianjie(real(X[1]),imag(X[1]),vp,vs,vf,denf,den,R,w,n,MTOOL,R_TOOL);
    f2=Funcxianjie(real(X[2]),imag(X[2]),vp,vs,vf,denf,den,R,w,n,MTOOL,R_TOOL);
    f3=Funcxianjie(real(X[3]),imag(X[3]),vp,vs,vf,denf,den,R,w,n,MTOOL,R_TOOL);
    f4=Funcxianjie(real(X[4]),imag(X[4]),vp,vs,vf,denf,den,R,w,n,MTOOL,R_TOOL);
    inter[0]=(complex<float>(1.0/12.0,0)/h)*complex<float>(-25,0)*f0;
    Inter=(complex<float>(1.0/12.0,0)/h)*complex<float>(48,0);
    inter[1]=inter[0]+Inter*f1;
    inter[2]=inter[1]+(complex<float>(1.0/12.0,0)/h)*complex<float>(-36,0)*f2;
    inter[3]=inter[2]+(complex<float>(1.0/12.0,0)/h)*complex<float>(16,0)*f3;
    inter[4]=inter[3]+(complex<float>(1.0/12.0,0)/h)*complex<float>(-3,0)*f4;

    m=inter[4];
    return m;
}

complex<float> CalcMethod::interposeW(float x,float y ,float vp,float vs,float vf,float denf,float den,float R,float w,int n,float MTOOL,float R_TOOL,float df)
{
    complex<float> m,inter[8],Inter;
    float W[5];
    int i,dff=df/5.;
    for(i=0;i<5;i++)
    {
        W[i]=w+(i)*dff;
    }
    complex<float> f0,f1,f2,f3,f4;
    f0=Funcxianjie(x,y,vp,vs,vf,denf,den,R,W[0],n,MTOOL,R_TOOL);
    f1=Funcxianjie(x,y,vp,vs,vf,denf,den,R,W[1],n,MTOOL,R_TOOL);
    f2=Funcxianjie(x,y,vp,vs,vf,denf,den,R,W[2],n,MTOOL,R_TOOL);
    f3=Funcxianjie(x,y,vp,vs,vf,denf,den,R,W[3],n,MTOOL,R_TOOL);
    f4=Funcxianjie(x,y,vp,vs,vf,denf,den,R,W[4],n,MTOOL,R_TOOL);

    inter[0]=complex<float>(1.0/(12.0*dff),0)*complex<float>(-25,0)*f0;
    Inter=complex<float>(48,0)*f1;
    inter[1]=inter[0]+complex<float>(1.0/(12.0*dff),0)*complex<float>(48,0)*f1;
    inter[2]=inter[1]+complex<float>(1.0/(12.0*dff),0)*complex<float>(-36,0)*f2;
    inter[3]=inter[2]+complex<float>(1.0/(12.0*dff),0)*complex<float>(16,0)*f3;
    inter[4]=inter[3]+complex<float>(1.0/(12.0*dff),0)*complex<float>(-3,0)*f4;
    m=inter[4];
    return m;
}
//
void CalcMethod::G_L(float r,float r_,float k,float w,float denf,float gl[2][2])
{
    float space=0.1524;
    gl[0][0]=r*r/(r_*r_)*cos(k*space);
    gl[0][1]=-r*r/(r_*r_)*sin(k*space)*k/(denf*w*w);
    gl[1][0]=denf*w*w*sin(k*space)/k;
    gl[1][1]=cos(k*space);

}
//求逆矩阵
void CalcMethod::MM_(float M[2][2],float M_[2][2])
{
    float x,x00,x01,x10,x11;
    x=M[0][0]*M[1][1]+M[0][1]*M[1][0];
    x00=M[1][1]/x;
    x01=M[1][0]/x;
    x10=M[0][1]/x;
    x11=M[0][0]/x;
    M_[0][0]=x00;
    M_[0][1]=x01;
    M_[1][0]=x10;
    M_[1][1]=x11;

}
//矩阵乘
void CalcMethod::M_M(float M1[2][2],float M2[2][2],float M[2][2])
{
    float x00,x01,x10,x11;
    x00=M1[0][0]*M2[0][0]+M1[0][1]*M2[1][0];
    x01=M1[0][0]*M2[0][1]+M1[0][1]*M2[1][1];
    x10=M1[1][0]*M2[0][0]+M1[1][1]*M2[1][0];
    x11=M1[1][0]*M2[0][1]+M1[1][1]*M2[1][1];
    M[0][0]=x00;
    M[0][1]=x01;
    M[1][0]=x10;
    M[1][1]=x11;
}
//
complex<float> CalcMethod::initialvalue(float x,float y ,float vp,float vs,float vf,float denf,float den,float R,float w,int n,float MTOOL,float R_TOOL,float df)
{
    complex<float> k;
    float W;
    W=df;//w/10;
    k=complex<float>(x,y);
    complex<float> difK,difW;
    difK=interposek(x,y ,vp,vs,vf,denf,den,R,w,n,MTOOL,R_TOOL);
    difW=interposeW(x,y ,vp,vs,vf,denf,den,R,w,n,MTOOL,R_TOOL,W);
    complex<float>temp=complex<float>(W,0)*difW/difK;
    k=k-temp;
    return k;
}