frechet.c

su 的源代码库

   bm.r = auxm3 * cos(angle);
   bm.i = auxm3 * sin(angle);

   /* bm * I */
   bmI.r = -bm.i;
   bmI.i = bm.r;
   
   /* computing phase-shift matrix */
   wThick.r = wC.r * (-2 * thick[0]);
   wThick.i = wC.i * (-2 * thick[0]);

   /* cexp (amI * wThick) */
   auxm1 = amI.r * wThick.r - amI.i * wThick.i;
   auxm2 = amI.r * wThick.i + amI.i * wThick.r;
   E[0][0].r = exp(auxm1) * cos(auxm2);
   E[0][0].i = exp(auxm1) * sin(auxm2);

   /* cexp((amI + bmI) * (wThick * .5)) */
   auxm1 = amI.r + bmI.r;
   auxm2 = amI.i + bmI.i;
   auxm3 = .5 * (auxm1 * wThick.r - auxm2 * wThick.i);
   auxm4 = .5 * (auxm1 * wThick.i + auxm2 * wThick.r);
   E[0][1].r = exp(auxm3) * cos(auxm4);
   E[0][1].i = exp(auxm3) * sin(auxm4);
   
   E[1][0] = E[0][1];

   /* cexp (bmI * wThick) */
   auxm1 = bmI.r * wThick.r - bmI.i * wThick.i;
   auxm2 = bmI.r * wThick.i + bmI.i * wThick.r;
   E[1][1].r = exp(auxm1) * cos(auxm2);
   E[1][1].i = exp(auxm1) * sin(auxm2);

   /* applying phase-shift */
   mT[0][0].r = mB[0][0].r * E[0][0].r - mB[0][0].i * E[0][0].i;
   mT[0][0].i = mB[0][0].r * E[0][0].i + mB[0][0].i * E[0][0].r;
   mT[0][1].r = mB[0][1].r * E[0][1].r - mB[0][1].i * E[0][1].i;
   mT[0][1].i = mB[0][1].r * E[0][1].i + mB[0][1].i * E[0][1].r;
   mT[1][0].r = mB[1][0].r * E[1][0].r - mB[1][0].i * E[1][0].i;
   mT[1][0].i = mB[1][0].r * E[1][0].i + mB[1][0].i * E[1][0].r;
   mT[1][1].r = mB[1][1].r * E[1][1].r - mB[1][1].i * E[1][1].i;
   mT[1][1].i = mB[1][1].r * E[1][1].i + mB[1][1].i * E[1][1].r;

   /* copying to matrix rm */
   rm[0][0] = mT[0][0]; rm[0][1] = mT[0][1];
   rm[1][0] = mT[1][0]; rm[1][1] = mT[1][1];

   /* applying phase-shift in the FRECHET derivatives as well */
   for (iDer = 0; iDer < numberPar; iDer++)
   {
      /* Vp -> Vs -> rho, if all active */
      for (i = 0, iL = MIN(lim[0], 2); i < limRange; iL++, i++) 
      {
	 aux = DmB[0][limRange * iDer + iL][0].r * E[0][0].r - 
	       DmB[0][limRange * iDer + iL][0].i * E[0][0].i;
	 DmB[0][limRange * iDer + iL][0].i = 
	       DmB[0][limRange * iDer + iL][0].r * E[0][0].i + 
	       DmB[0][limRange * iDer + iL][0].i * E[0][0].r;
	 DmB[0][limRange * iDer + iL][0].r = aux;
	 aux = DmB[0][limRange * iDer + iL][1].r * E[0][1].r - 
	       DmB[0][limRange * iDer + iL][1].i * E[0][1].i;
	 DmB[0][limRange * iDer + iL][1].i = 
	       DmB[0][limRange * iDer + iL][1].r * E[0][1].i + 
	       DmB[0][limRange * iDer + iL][1].i * E[0][1].r;
	 DmB[0][limRange * iDer + iL][1].r = aux;
	 aux = DmB[0][limRange * iDer + iL][2].r * E[1][0].r - 
	       DmB[0][limRange * iDer + iL][2].i * E[1][0].i;
	 DmB[0][limRange * iDer + iL][2].i = 
	       DmB[0][limRange * iDer + iL][2].r * E[1][0].i + 
	       DmB[0][limRange * iDer + iL][2].i * E[1][0].r;
	 DmB[0][limRange * iDer + iL][2].r = aux;
	 aux = DmB[0][limRange * iDer + iL][3].r * E[1][1].r - 
	       DmB[0][limRange * iDer + iL][3].i * E[1][1].i;
	 DmB[0][limRange * iDer + iL][3].i = 
	       DmB[0][limRange * iDer + iL][3].r * E[1][1].i + 
	       DmB[0][limRange * iDer + iL][3].i * E[1][1].r;
	 DmB[0][limRange * iDer + iL][3].r = aux;
      }
   }
}
/*                                                              */
/*  Function frechetRm()                                        */
/*                                                              */
/*  Computing Frechet derivatives for the reflectivity          */
/*  matrix rm                                                   */
/*                                                              */
/*  Input parameters:                                           */
/*  E[2][2]................phase shift matrix                   */
/*  tD[2][2]...............downgoing transmission coefficients  */
/*  rU[2][2]...............upgoing reflection coefficients      */
/*  mT.....................partial reflectivity matrix          */
/*  inv....................(I - mT * rU)^-1                     */
/*  wThick.................wC * thickness                       */
/*  am.....................p-wave vertical slowness             */
/*  bm.....................s-wave vertical slowness             */
/*  iL.....................points to the current layer          */
/*                                                              */
/*  Output parameters:                                          */
/*  DmB....................Frechet derivatives of the           */
/*                         reflectivity matrix at level iL      */
/*                                                              */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
void frechetRm(complex E[2][2], complex tD[2][2], complex tUinv[2][2], 
	       complex mTtD[2][2], complex mT[2][2], complex rU[2][2], 
	       complex inv[2][2], complex wThick, complex am, 
	       complex bm, int iL)
{
   /* declaration of variables */
   int iDer, k, kDer=0, kPnt=0;           /* counters */
   complex sInv;                          /* 1 / Slowness */
   static complex DrD[2][2], DtD[2][2];   /* coefficient derivatives */
   static complex DrU[2][2], DtU[2][2];   /* coefficient derivatives */
   static complex factor1[2][2], factor2[2][2], factor3[2][2];
   static complex A[2][2], B[2][2], C[2][2];
                                          /* auxiliar variables */
   static complex invmTtD[2][2];          /* inv * mT * tD */
   
   /* inv * mTtD */
   matMult(invmTtD, inv, mTtD);

   /* pointer to previous derivatives */
   for (k = MAX(iL - 1, lim[0]), kDer = 2; k < lim[1]; k++)
   {
      /* considering special cases */
      if (k == iL - 1)
      {
	 /* differentiating reflection and transmition coefficients */
	 frechetRTd(iL - 1, iL, iL - 1);
	 frechetRTu(iL - 1, iL, iL - 1);
	 
	 for (iDer = 0; iDer < numberPar; iDer++)
	 {
	    /* Vp -> Vs -> rho, if all active */
	    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);
	    /* mT * DrU */
	    matMult(A, mT, DrU);
	    /* A * mTtD */
	    matMult(B, A, invmTtD);
	    /* C * tUinv */
	    matMult(factor2, tUinv, B);
	    
	    /* mT * DtD */
	    matMult(A, mT, DtD);
	    /* tUinv * A */
	    matMult(factor3, tUinv, A);
	    
	    /* compounding factors */
	    /*      [ 0  1 ]       */
	    /*      [ 2  3 ]       */
	    DmB[iL - 1][limRange * iDer + 0][0].r = DrD[0][0].r + 
	       factor1[0][0].r + factor2[0][0].r + factor3[0][0].r;
	    DmB[iL - 1][limRange * iDer + 0][0].i = DrD[0][0].i + 
	       factor1[0][0].i + factor2[0][0].i + factor3[0][0].i;
	    DmB[iL - 1][limRange * iDer + 0][1].r = DrD[0][1].r + 
	       factor1[0][1].r + factor2[0][1].r + factor3[0][1].r;
	    DmB[iL - 1][limRange * iDer + 0][1].i = DrD[0][1].i + 
	       factor1[0][1].i + factor2[0][1].i + factor3[0][1].i;
	    DmB[iL - 1][limRange * iDer + 0][2].r = DrD[1][0].r + 
	       factor1[1][0].r + factor2[1][0].r + factor3[1][0].r;
	    DmB[iL - 1][limRange * iDer + 0][2].i = DrD[1][0].i + 
	       factor1[1][0].i + factor2[1][0].i + factor3[1][0].i;
	    DmB[iL - 1][limRange * iDer + 0][3].r = DrD[1][1].r + 
	       factor1[1][1].r + factor2[1][1].r + factor3[1][1].r;
	    DmB[iL - 1][limRange * iDer + 0][3].i = DrD[1][1].i + 
	       factor1[1][1].i +  factor2[1][1].i + factor3[1][1].i;

	 }
      }
      else if (k == iL)
      {
	 /* differentiating reflection and transmition coefficients */
	 frechetRTd(iL - 1, iL, iL);
	 frechetRTu(iL - 1, iL, iL);
	 
	 /* 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;