frechetslave.c

su 的源代码库

	 /* saving flags */
	 vpF = vpFrechet;
	 vsF = vsFrechet;
	 rhoF = rhoFrechet;
	 for (iDer = 0; iDer < numberPar; iDer++)
	 {
	    /* Vp -> Vs -> rho */
	    DrD[0][0] = coeffDFr[iDer][0]; DrD[0][1] = coeffDFr[iDer][1]; 
	    DrD[1][0] = coeffDFr[iDer][2]; DrD[1][1] = coeffDFr[iDer][3];
	 
	    DtD[0][0] = coeffDFr[iDer][4]; DtD[0][1] = coeffDFr[iDer][5]; 
	    DtD[1][0] = coeffDFr[iDer][6]; DtD[1][1] = coeffDFr[iDer][7];
	    
	    DrU[0][0] = coeffUFr[iDer][0]; DrU[0][1] = coeffUFr[iDer][1]; 
	    DrU[1][0] = coeffUFr[iDer][2]; DrU[1][1] = coeffUFr[iDer][3];
	    
	    DtU[0][0] = coeffUFr[iDer][4]; DtU[0][1] = coeffUFr[iDer][5]; 
	    DtU[1][0] = coeffUFr[iDer][6]; DtU[1][1] = coeffUFr[iDer][7];
	    
	    /* DtU * invmTtD */
	    matMult(factor1, DtU, invmTtD);
	    
	    /* applying phase-shift on D (MB_iL+1) / D (iL+1) */
	    /*      [ 0  1 ]       */
	    /*      [ 2  3 ]       */
	    A[0][0].r = DmB[iL][limRange * iDer + 0][0].r * E[0][0].r
  	              - DmB[iL][limRange * iDer + 0][0].i * E[0][0].i;
	    A[0][0].i = DmB[iL][limRange * iDer + 0][0].r * E[0][0].i 
	              + DmB[iL][limRange * iDer + 0][0].i * E[0][0].r;
	    A[0][1].r = DmB[iL][limRange * iDer + 0][1].r * E[0][1].r 
	              - DmB[iL][limRange * iDer + 0][1].i * E[0][1].i;
	    A[0][1].i = DmB[iL][limRange * iDer + 0][1].r * E[0][1].i 
	              + DmB[iL][limRange * iDer + 0][1].i * E[0][1].r;
	    A[1][0].r = DmB[iL][limRange * iDer + 0][2].r * E[1][0].r 
	              - DmB[iL][limRange * iDer + 0][2].i * E[1][0].i;
	    A[1][0].i = DmB[iL][limRange * iDer + 0][2].r * E[1][0].i 
	              + DmB[iL][limRange * iDer + 0][2].i * E[1][0].r;
	    A[1][1].r = DmB[iL][limRange * iDer + 0][3].r * E[1][1].r 
	              - DmB[iL][limRange * iDer + 0][3].i * E[1][1].i;
	    A[1][1].i = DmB[iL][limRange * iDer + 0][3].r * E[1][1].i 
	              + DmB[iL][limRange * iDer + 0][3].i * E[1][1].r;
	    
	    /* things now are different if the derivative is being */
	    /* taken with respect to P_wave, S_wave velocities or */
	    /* densities */
	    if (vpFrechet) /* P-wave */
	    {
	       /* mT * [i wThick] */
	       /* note that wThick includes the - sign */
	       auxm1 = -wThick.i; auxm2 = wThick.r;
	       B[0][0].r = mT[0][0].r * auxm1 - mT[0][0].i * auxm2;
	       B[0][0].i = mT[0][0].r * auxm2 + mT[0][0].i * auxm1;
	       B[0][1].r = 0.5 * (mT[0][1].r * auxm1 - mT[0][1].i * auxm2);
	       B[0][1].i = 0.5 * (mT[0][1].r * auxm2 + mT[0][1].i * auxm1);
	       B[1][0].r = 0.5 * (mT[1][0].r * auxm1 - mT[1][0].i * auxm2);
	       B[1][0].i = 0.5 * (mT[1][0].r * auxm2 + mT[1][0].i * auxm1);
	       B[1][1].r = mT[1][1].r * auxm1 - mT[1][1].i * auxm2;
	       B[1][1].i = mT[1][1].r * auxm2 + mT[1][1].i * auxm1;
	       
	       /* 1 / am */
	       aux = am.r * am.r + am.i * am.i;
	       sInv.r = am.r / aux;
	       sInv.i = -am.i / aux;
	       
	       /* 1 / am * derFactor[0] */
	       auxm1 = sInv.r * derFactor[iL][0].r - 
		       sInv.i * derFactor[iL][0].i;
	       auxm2 = sInv.r * derFactor[iL][0].i + 
		       sInv.i * derFactor[iL][0].r;
	       
	       C[0][0].r = A[0][0].r - (B[0][0].r * auxm1 - B[0][0].i * auxm2);
	       C[0][0].i = A[0][0].i - (B[0][0].r * auxm2 + B[0][0].i * auxm1);
	       C[0][1].r = A[0][1].r - (B[0][1].r * auxm1 - B[0][1].i * auxm2);
	       C[0][1].i = A[0][1].i - (B[0][1].r * auxm2 + B[0][1].i * auxm1);
	       C[1][0].r = A[1][0].r - (B[1][0].r * auxm1 - B[1][0].i * auxm2);
	       C[1][0].i = A[1][0].i - (B[1][0].r * auxm2 + B[1][0].i * auxm1);
	       C[1][1] = A[1][1];
	       vpFrechet = 0;
	    }
	    else if (vsFrechet) /* S-wave */
	    {
	       /* mT * [i -wThick] */
	       /* note that wThick includes the - sign */
	       auxm1 = -wThick.i; auxm2 = wThick.r;
	       B[0][0].r = mT[0][0].r * auxm1 - mT[0][0].i * auxm2;
	       B[0][0].i = mT[0][0].r * auxm2 + mT[0][0].i * auxm1;
	       B[0][1].r = 0.5 * (mT[0][1].r * auxm1 - mT[0][1].i * auxm2);
	       B[0][1].i = 0.5 * (mT[0][1].r * auxm2 + mT[0][1].i * auxm1);
	       B[1][0].r = 0.5 * (mT[1][0].r * auxm1 - mT[1][0].i * auxm2);
	       B[1][0].i = 0.5 * (mT[1][0].r * auxm2 + mT[1][0].i * auxm1);
	       B[1][1].r = mT[1][1].r * auxm1 - mT[1][1].i * auxm2;
	       B[1][1].i = mT[1][1].r * auxm2 + mT[1][1].i * auxm1;

	       /* 1 / bm */
	       aux = bm.r * bm.r + bm.i * bm.i;
	       sInv.r = bm.r / aux;
	       sInv.i = -bm.i / aux;
	       
	       /* 1 / Slowness * derFactor[1] */
	       auxm1 = sInv.r * derFactor[iL][1].r - 
		       sInv.i * derFactor[iL][1].i;
	       auxm2 = sInv.r * derFactor[iL][1].i + 
		       sInv.i * derFactor[iL][1].r;
	    
	       C[0][0] = A[0][0];
	       C[0][1].r = A[0][1].r - (B[0][1].r * auxm1 - B[0][1].i * auxm2);
	       C[0][1].i = A[0][1].i - (B[0][1].r * auxm2 + B[0][1].i * auxm1);
	       C[1][0].r = A[1][0].r - (B[1][0].r * auxm1 - B[1][0].i * auxm2);
	       C[1][0].i = A[1][0].i - (B[1][0].r * auxm2 + B[1][0].i * auxm1);
	       C[1][1].r = A[1][1].r - (B[1][1].r * auxm1 - B[1][1].i * auxm2);
	       C[1][1].i = A[1][1].i - (B[1][1].r * auxm2 + B[1][1].i * auxm1);
	       vsFrechet = 0;
	    }
	    else if (rhoFrechet) /* Density */
	    {
	       C[0][0] = A[0][0];
	       C[0][1] = A[0][1];
	       C[1][0] = A[1][0];
	       C[1][1] = A[1][1];
	       rhoFrechet = 0;
	    }
	 
	    /* C * rU */
	    matMult(A, C, rU);
	    /* mT * DrU */
	    matMult(B, mT, DrU);
	    A[0][0].r += B[0][0].r;
	    A[0][0].i += B[0][0].i;
	    A[0][1].r += B[0][1].r;
	    A[0][1].i += B[0][1].i;
	    A[1][0].r += B[1][0].r;
	    A[1][0].i += B[1][0].i;
	    A[1][1].r += B[1][1].r;
	    A[1][1].i += B[1][1].i;
	    
	    /* A * invmTtD */
	    matMult(B, A, invmTtD);
	    /* tUinv * B */
	    matMult(factor2, tUinv, B);
	    
	    /* C * tD */
	    matMult(A, C, tD);
	    /* mT * DtD */
	    matMult(B, mT, DtD);
	    
	    A[0][0].r += B[0][0].r; A[0][0].i += B[0][0].i;
	    A[1][0].r += B[1][0].r; A[1][0].i += B[1][0].i;
	    A[0][1].r += B[0][1].r; A[0][1].i += B[0][1].i;
	    A[1][1].r += B[1][1].r; A[1][1].i += B[1][1].i;
	    
	    /* tUinv * A */
	    matMult(factor3, tUinv, A);
	    
	    /* compounding factors */
	    /*      [ 0  1 ]       */
	    /*      [ 2  3 ]       */
	    DmB[iL - 1][limRange * iDer + 1][0].r = DrD[0][0].r + 
	       factor1[0][0].r + factor2[0][0].r + factor3[0][0].r;
	    DmB[iL - 1][limRange * iDer + 1][0].i = DrD[0][0].i + 
	       factor1[0][0].i + factor2[0][0].i + factor3[0][0].i;
	    DmB[iL - 1][limRange * iDer + 1][1].r = DrD[0][1].r + 
	       factor1[0][1].r + factor2[0][1].r + factor3[0][1].r;
	    DmB[iL - 1][limRange * iDer + 1][1].i = DrD[0][1].i + 
	       factor1[0][1].i + factor2[0][1].i + factor3[0][1].i;
	    DmB[iL - 1][limRange * iDer + 1][2].r = DrD[1][0].r + 
	       factor1[1][0].r + factor2[1][0].r + factor3[1][0].r;
	    DmB[iL - 1][limRange * iDer + 1][2].i = DrD[1][0].i + 
	       factor1[1][0].i + factor2[1][0].i + factor3[1][0].i;
	    DmB[iL - 1][limRange * iDer + 1][3].r = DrD[1][1].r + 
	       factor1[1][1].r + factor2[1][1].r + factor3[1][1].r;
	    DmB[iL - 1][limRange * iDer + 1][3].i = DrD[1][1].i + 
	       factor1[1][1].i + factor2[1][1].i + factor3[1][1].i;
	 }
	 /* restoring */
	 vpFrechet = vpF;
	 vsFrechet = vsF;
	 rhoFrechet = rhoF;
      }
      else
      {
	 /* applying phase-shift on D (MB_j+1) / D (j...on) */
	 /*      [ 0  1 ]       */
	 /*      [ 2  3 ]       */
	 for (iDer = 0; iDer < numberPar; iDer++)
	 {
	    /* Vp -> Vs -> rho, if all active */
	    kPnt = ((lim[0] - iL) >= (2) ? (kDer) : (kDer - 1));
	    A[0][0].r = DmB[iL][limRange * iDer + kPnt][0].r * E[0][0].r
	              - DmB[iL][limRange * iDer + kPnt][0].i * E[0][0].i;
	    A[0][0].i = DmB[iL][limRange * iDer + kPnt][0].r * E[0][0].i 
	              + DmB[iL][limRange * iDer + kPnt][0].i * E[0][0].r;
	    A[0][1].r = DmB[iL][limRange * iDer + kPnt][1].r * E[0][1].r 
	              - DmB[iL][limRange * iDer + kPnt][1].i * E[0][1].i;
	    A[0][1].i = DmB[iL][limRange * iDer + kPnt][1].r * E[0][1].i 
	              + DmB[iL][limRange * iDer + kPnt][1].i * E[0][1].r;
	    A[1][0].r = DmB[iL][limRange * iDer + kPnt][2].r * E[1][0].r 
	              - DmB[iL][limRange * iDer + kPnt][2].i * E[1][0].i;
	    A[1][0].i = DmB[iL][limRange * iDer + kPnt][2].r * E[1][0].i 
	              + DmB[iL][limRange * iDer + kPnt][2].i * E[1][0].r;
	    A[1][1].r = DmB[iL][limRange * iDer + kPnt][3].r * E[1][1].r 
	              - DmB[iL][limRange * iDer + kPnt][3].i * E[1][1].i;
	    A[1][1].i = DmB[iL][limRange * iDer + kPnt][3].r * E[1][1].i 
	              + DmB[iL][limRange * iDer + kPnt][3].i * E[1][1].r;
 
	    /* tUinv * A */
	    matMult(C, tUinv, A);
	    
	    /* rU * invmTtD */
	    matMult(B, rU, invmTtD);
	    B[0][0].r += tD[0][0].r; B[0][0].i += tD[0][0].i;
	    B[0][1].r += tD[0][1].r; B[0][1].i += tD[0][1].i;
	    B[1][0].r += tD[1][0].r; B[1][0].i += tD[1][0].i;
	    B[1][1].r += tD[1][1].r; B[1][1].i += tD[1][1].i;
	    
	    /* C * B */
	    matMult(factor1, C, B);
	    
	    /* compounding factors */
	    /*      [ 0  1 ]       */
	    /*      [ 2  3 ]       */
	    DmB[iL - 1][limRange * iDer + kDer][0].r = factor1[0][0].r;
	    DmB[iL - 1][limRange * iDer + kDer][0].i = factor1[0][0].i;
	    DmB[iL - 1][limRange * iDer + kDer][1].r = factor1[0][1].r;
	    DmB[iL - 1][limRange * iDer + kDer][1].i = factor1[0][1].i;
	    DmB[iL - 1][limRange * iDer + kDer][2].r = factor1[1][0].r;
	    DmB[iL - 1][limRange * iDer + kDer][2].i = factor1[1][0].i;
	    DmB[iL - 1][limRange * iDer + kDer][3].r = factor1[1][1].r;
	    DmB[iL - 1][limRange * iDer + kDer][3].i = factor1[1][1].i;
	 }
	 kDer ++;
      }
   }
}
/*                                                              */
/*  Function matMult()                                          */
/*                                                              */
/*  Multiplication of two complex 2x2 matrices                  */
/*                                                              */
/*  Input parameters:                                           */
/*  A[2][2]................input matrix                         */
/*  B[2][2]................input matrix                         */
/*                                                              */
/*  Output parameters:                                          */
/*  C[2][2]................A * B                                */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
void matMult(complex A[2][2], complex B[2][2], complex C[2][2])
{
      A[0][0].r = B[0][0].r * C[0][0].r - B[0][0].i * C[0][0].i
  	        + B[0][1].r * C[1][0].r - B[0][1].i * C[1][0].i;
      A[0][0].i = B[0][0].r * C[0][0].i + B[0][0].i * C[0][0].r
	        + B[0][1].r * C[1][0].i + B[0][1].i * C[1][0].r;
      A[0][1].r = B[0][0].r * C[0][1].r - B[0][0].i * C[0][1].i
	        + B[0][1].r * C[1][1].r - B[0][1].i * C[1][1].i;
      A[0][1].i = B[0][0].r * C[0][1].i + B[0][0].i * C[0][1].r
	        + B[0][1].r * C[1][1].i + B[0][1].i * C[1][1].r;
      A[1][0].r = B[1][0].r * C[0][0].r - B[1][0].i * C[0][0].i
                + B[1][1].r * C[1][0].r - B[1][1].i * C[1][0].i;
      A[1][0].i = B[1][0].r * C[0][0].i + B[1][0].i * C[0][0].r
      	        + B[1][1].r * C[1][0].i + B[1][1].i * C[1][0].r;
      A[1][1].r = B[1][0].r * C[0][1].r - B[1][0].i * C[0][1].i
      	        + B[1][1].r * C[1][1].r - B[1][1].i * C[1][1].i;
      A[1][1].i = B[1][0].r * C[0][1].i + B[1][0].i * C[0][1].r
      	        + B[1][1].r * C[1][1].i + B[1][1].i * C[1][1].r;
}
/*                                                              */
/*  Function frechetRTD()                                       */
/*                                                              */
/*  Numerical computation of the derivatives of the             */
/*  reflection and transmission coefficients of downgoing       */
/*  incident waves with respect to the model parameters         */
/*                                                              */
/*  Input parameters:                                           */
/*  i1.....................points to the first layer            */
/*  i2.....................points to the second layer           */
/*  iF.....................indicates if the derivatives are with*/
/*                         respect to the first or second layer */
/*                                                              */
/*  Output parameters:                                          */
/*  coeffDFr...............Derivatives of the elastic           */
/*                         reflection coefficients with respect */
/*                         to Vp, Vs and rho                    */
/*                                                              */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
void frechetRTd(int i1, int i2, int iF)
{ 
   /* declaration of variables */
   int iDer;                       /* counter */
   float rho1=0, rho2=0;           /* densities */
   float rho1rho2;                 /* rho1 * rho2 */
   float DIV=0;                    /* derivative denominator */
   complex aux1, aux2, aux3, c;    
   complex cuu, caux3, aux1aux2;   /* auxiliar quantities */
   complex a1, a2, b1, b2;         /* all vertical slownesses */
   complex a1b1, a1b2, a2b1, a2b2; /* auxiliar quantities */
   complex a1a2b1b2;               /* auxiliar quantities */
   complex d1d, d2d;               /* auxiliar quantities */
   complex dpda, dpdb;             /* auxiliar quantities */
   complex PSlow1, PSlow2;         /* P-wave slownesses */
   complex SSlow1, SSlow2;         /* S-wave slownesses */
   complex S2Vel1, S2Vel2;         /* S-wave velocities */

   /* iDer indicates what parameter is active */
   /* bookkeeping the parameters */
   vpF = vpFrechet;
   vsF = vsFrechet;
   rhoF = rhoFrechet;

   for (iDer = 0; iDer < numberPar; iDer++)
   {
      /* Vp -> Vs -> rho, if all active */
      if (vpFrechet)
      {
	 SSlow1 = SSlowness[i1][0];
	 SSlow2 = SSlowness[i2][0];
	 S2Vel1 = S2Velocity[i1][0];
	 S2Vel2 = S2Velocity[i2][0];
	 rho1 = rho[i1];
	 rho2 = rho[i2];
	 if (iF == i1)
	 {
	    PSlow1 = PSlowness[