frechet.c

su 的源代码库

	    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[i1][1];
	    PSlow2 = PSlowness[i2][0];
	    DIV = (FACTOR - 1.) * alpha[i1];
	 }
	 else if (iF == i2)
	 {
	    PSlow1 = PSlowness[i1][0];
	    PSlow2 = PSlowness[i2][1];
	    DIV = (FACTOR - 1.) * alpha[i2];
	 }
	 vpFrechet = 0;
      }
      else if (vsFrechet)
      {
	 PSlow1 = PSlowness[i1][0];
	 PSlow2 = PSlowness[i2][0];
	 rho1 = rho[i1];
	 rho2 = rho[i2];
	 if (iF == i1)
	 {
	    SSlow1 = SSlowness[i1][1];
	    SSlow2 = SSlowness[i2][0];
	    S2Vel1 = S2Velocity[i1][1];
	    S2Vel2 = S2Velocity[i2][0];
	    DIV = (FACTOR - 1.) * beta[i1];
	 }
	 else if (iF == i2)
	 {
	    SSlow1 = SSlowness[i1][0];
	    SSlow2 = SSlowness[i2][1];
	    S2Vel1 = S2Velocity[i1][0];
	    S2Vel2 = S2Velocity[i2][1];
	    DIV = (FACTOR - 1.) * beta[i2];
	 }
	 vsFrechet = 0;
      }
      else if (rhoFrechet)
      {
	 PSlow1 = PSlowness[i1][0];
	 PSlow2 = PSlowness[i2][0];
	 SSlow1 = SSlowness[i1][0];
	 SSlow2 = SSlowness[i2][0];
	 S2Vel1 = S2Velocity[i1][0];
	 S2Vel2 = S2Velocity[i2][0];
	 if (iF == i1)
	 {
	    rho1 = FACTOR * rho[i1];
	    rho2 = rho[i2];
	    DIV = (FACTOR - 1.) * rho[i1];
	 }
	 else if (iF == i2)
	 {
	    rho1 = rho[i1];
	    rho2 = FACTOR * rho[i2];
	    DIV = (FACTOR - 1.) * rho[i2];
	 }
	 rhoFrechet = 0;
      }
      
      /* square-root of PSlowness^2 - uuC */    
      auxm1 = PSlow1.r - uuC.r;
      auxm2 = PSlow1.i - uuC.i;
      auxm3 = sqrt(auxm1 * auxm1 + auxm2 * auxm2);
      auxm3 = sqrt(auxm3);
      angle = atan2(auxm2, auxm1) / 2;
      a1.r = auxm3 * cos(angle);
      a1.i = auxm3 * sin(angle);
      
      /* square-root of SSlowness^2 - uuC */    
      auxm1 = SSlow1.r - uuC.r;
      auxm2 = SSlow1.i - uuC.i;
      auxm3 = sqrt(auxm1 * auxm1 + auxm2 * auxm2);
      auxm3 = sqrt(auxm3);
      angle = atan2(auxm2, auxm1) / 2;
      b1.r = auxm3 * cos(angle);
      b1.i = auxm3 * sin(angle);
      
      /* square-root of PSlowness^2 - uuC */    
      auxm1 = PSlow2.r - uuC.r;
      auxm2 = PSlow2.i - uuC.i;
      auxm3 = sqrt(auxm1 * auxm1 + auxm2 * auxm2);
      auxm3 = sqrt(auxm3);
      angle = atan2(auxm2, auxm1) / 2;
      a2.r = auxm3 * cos(angle);
      a2.i = auxm3 * sin(angle);
      
      /* square-root of SSlowness^2 - uuC */    
      auxm1 = SSlow2.r - uuC.r;
      auxm2 = SSlow2.i - uuC.i;
      auxm3 = sqrt(auxm1 * auxm1 + auxm2 * auxm2);
      auxm3 = sqrt(auxm3);
      angle = atan2(auxm2, auxm1) / 2;
      b2.r = auxm3 * cos(angle);
      b2.i = auxm3 * sin(angle);
      
      rho1rho2 = rho1 * rho2;
      c.r = 2 * (rho1 * S2Vel1.r - rho2 * S2Vel2.r);
      c.i = 2 * (rho1 * S2Vel1.i - rho2 * S2Vel2.i);
      
      cuu.r = c.r * uuC.r - c.i * uuC.i;
      cuu.i = c.r * uuC.i + c.i * uuC.r;
      
      aux1.r = cuu.r - (rho1 - rho2);
      aux1.i = cuu.i;
      
      aux2.r = cuu.r + rho2;
      aux2.i = cuu.i;
      
      aux1aux2.r = aux1.r * aux2.r - aux1.i * aux2.i;
      aux1aux2.i = aux1.r * aux2.i + aux1.i * aux2.r;
      
      aux3.r = cuu.r - rho1;
      aux3.i = cuu.i;
      
      caux3.r = c.r * aux3.r - c.i * aux3.i;
      caux3.i = c.r * aux3.i + c.i * aux3.r;
      
      a1b1.r = a1.r * b1.r - a1.i * b1.i;
      a1b1.i = a1.r * b1.i + a1.i * b1.r;
      a1b2.r = a1.r * b2.r - a1.i * b2.i;
      a1b2.i = a1.r * b2.i + a1.i * b2.r;
      a2b1.r = a2.r * b1.r - a2.i * b1.i;
      a2b1.i = a2.r * b1.i + a2.i * b1.r;
      a2b2.r = a2.r * b2.r - a2.i * b2.i;