frechet.c

su 的源代码库

      b2.r = auxm3 * cos(angle);
      b2.i = auxm3 * sin(angle);
      
      /* computing auxiliary quantities */
      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;
      
      aux3.r = cuu.r - rho1;
      aux3.i = cuu.i;
      
      caux2.r = c.r * aux2.r - c.i * aux2.i;
      caux2.i = c.r * aux2.i + c.i * aux2.r;
      
      aux1aux3.r = aux1.r * aux3.r - aux1.i * aux3.i;
      aux1aux3.i = aux1.r * aux3.i + aux1.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;
      a2b2.i = a2.r * b2.i + a2.i * b2.r;
      a1a2b1b2.r = a1b1.r * a2b2.r - a1b1.i * a2b2.i;
      a1a2b1b2.i = a1b1.r * a2b2.i + a1b1.i * a2b2.r;
      
      /* computing D factors */
      auxm1 = aux2.r * aux2.r - aux2.i * aux2.i;
      auxm2 = 2 * aux2.r * aux2.i;
      auxm3 = a1b1.r * auxm1 - a1b1.i * auxm2;
      auxm4 = a1b1.r * auxm2 + a1b1.i * auxm1;
      d1u.r = auxm3 + a1b2.r * rho1rho2;
      d1u.i = auxm4 + a1b2.i * rho1rho2;
      
      auxm1 = aux1.r * aux1.r - aux1.i * aux1.i;
      auxm2 = 2 * aux1.r * aux1.i;
      auxm3 = auxm1 * uuC.r - auxm2 * uuC.i;
      auxm4 = auxm1 * uuC.i + auxm2 * uuC.r;
      d1u.r += auxm3;
      d1u.i += auxm4;
      
      auxm1 = aux3.r * aux3.r - aux3.i * aux3.i;
      auxm2 = 2 * aux3.r * aux3.i;
      auxm3 = a2b2.r * auxm1 - a2b2.i * auxm2;
      auxm4 = a2b2.r * auxm2 + a2b2.i * auxm1;
      d2u.r = auxm3 + a2b1.r * rho1rho2;
      d2u.i = auxm4 + a2b1.i * rho1rho2;
      
      auxm1 = c.r * cuu.r - c.i * cuu.i;
      auxm2 = c.r * cuu.i + c.i * cuu.r;
      auxm3 = a1a2b1b2.r * auxm1 - a1a2b1b2.i * auxm2;
      auxm4 = a1a2b1b2.r * auxm2 + a1a2b1b2.i * auxm1;
      d2u.r += auxm3;
      d2u.i += auxm4;
      
      /* more auxiliar quantities */
      /* d1u + d2u */
      dd.r = d1u.r + d2u.r;
      dd.i = d1u.i + d2u.i;
      
      /* 1 / dd */
      aux = dd.r * dd.r + dd.i * dd.i;
      dd.r = dd.r / aux;
      dd.i = -dd.i / aux;
      
      dpda.r = a2.r * dd.r - a2.i * dd.i;
      dpda.i = a2.r * dd.i + a2.i * dd.r;
      dpdb.r = b2.r * dd.r - b2.i * dd.i;
      dpdb.i = b2.r * dd.i + b2.i * dd.r;
      
      /* computing the coefficients - first reflection */
      auxm1 = d2u.r - d1u.r;
      auxm2 = d2u.i - d1u.i;
      /* (d2u - d1u) / (d1u + d2u) */
      coeffUFr[iDer][0].r = auxm1 * dd.r - auxm2 * dd.i;
      coeffUFr[iDer][0].i = auxm1 * dd.i + auxm2 * dd.r;    /* Rpp */
      
      /* computing the derivative */
      coeffUFr[iDer][0].r = (coeffUFr[iDer][0].r - coeffU[0].r) / DIV;
      coeffUFr[iDer][0].i = (coeffUFr[iDer][0].i - coeffU[0].i) / DIV;
      
      auxm1 = a1b1.r * caux2.r - a1b1.i * caux2.i;
      auxm2 = a1b1.r * caux2.i + a1b1.i * caux2.r;
      coeffUFr[iDer][1].r = aux1aux3.r + auxm1;
      coeffUFr[iDer][1].i = aux1aux3.i + auxm2;             /* Rsp */  
      
      coeffUFr[iDer][2].r = coeffUFr[iDer][1].r;
      coeffUFr[iDer][2].i = coeffUFr[iDer][1].i;            /* Rps */  
      
      auxm3 = -dpdb.r * uC2.r + dpdb.i * uC2.i;
      auxm4 = -dpdb.r * uC2.i - dpdb.i * uC2.r;
      
      aux = auxm3 * coeffUFr[iDer][1].r - auxm4 * coeffUFr[iDer][1].i;
      coeffUFr[iDer][1].i = auxm3 * coeffUFr[iDer][1].i + 
	                    auxm4 * coeffUFr[iDer][1].r;
      coeffUFr[iDer][1].r = aux;                            /* Rsp */
   
      /* computing the derivative */
      coeffUFr[iDer][1].r = (coeffUFr[iDer][1].r - coeffU[1].r) / DIV;
      coeffUFr[iDer][1].i = (coeffUFr[iDer][1].i - coeffU[1].i) / DIV;

      auxm3 = dpda.r * uC2.r - dpda.i * uC2.i;
      auxm4 = dpda.r * uC2.i + dpda.i * uC2.r;
      
      aux = auxm3 * coeffUFr[iDer][2].r - auxm4 * coeffUFr[iDer][2].i;
      coeffUFr[iDer][2].i = auxm3 * coeffUFr[iDer][2].i + 
	                    auxm4 * coeffUFr[iDer][2].r;
      coeffUFr[iDer][2].r = aux;                           /* Rps */
      
      /* computing the derivative */
      coeffUFr[iDer][2].r = (coeffUFr[iDer][2].r - coeffU[2].r) / DIV;
      coeffUFr[iDer][2].i = (coeffUFr[iDer][2].i - coeffU[2].i) / DIV;
      auxm1 = d2u.r - d1u.r;
      auxm2 = d2u.i - d1u.i;
      auxm3 = 2 * rho1rho2 * (a2b1.r - a1b2.r);
      auxm4 = 2 * rho1rho2 * (a2b1.i - a1b2.i);
      
      coeffUFr[iDer][3].r = auxm1 - auxm3;
      coeffUFr[iDer][3].i = auxm2 - auxm4;
      aux = coeffUFr[iDer][3].r * dd.r - coeffUFr[iDer][3].i * dd.i;
      coeffUFr[iDer][3].i = coeffUFr[iDer][3].r * dd.i + 
	                    coeffUFr[iDer][3].i * dd.r;
      coeffUFr[iDer][3].r = aux;                          /* Rss */
      
      /* computing the derivative */
      coeffUFr[iDer][3].r = (coeffUFr[iDer][3].r - coeffU[3].r) / DIV;
      coeffUFr[iDer][3].i = (coeffUFr[iDer][3].i - coeffU[3].i) / DIV;
      /* now transmition */
      auxm1 = b1.r * aux2.r - b1.i * aux2.i;
      auxm2 = b1.r * aux2.i + b1.i * aux2.r;
      auxm3 = b2.r * aux3.r - b2.i * aux3.i;
      auxm4 = b2.r * aux3.i + b2.i * aux3.r;
      coeffUFr[iDer][4].r = auxm1 - auxm3;     
      coeffUFr[iDer][4].i = auxm2 - auxm4;
      aux = 2 * rho2 * (dpda.r * coeffUFr[iDer][4].r - 
			dpda.i * coeffUFr[iDer][4].i);
      coeffUFr[iDer][4].i = 2 * rho2 * (dpda.r * coeffUFr[iDer][4].i + 
					dpda.i * coeffUFr[iDer][4].r);
      coeffUFr[iDer][4].r = aux;                            /* Tpp */
      
      /* computing the derivative */
      coeffUFr[iDer][4].r = (coeffUFr[iDer][4].r - coeffU[4].r) / DIV;
      coeffUFr[iDer][4].i = (coeffUFr[iDer][4].i - coeffU[4].i) / DIV;
      
      coeffUFr[iDer][5].r = aux1.r + a2b1.r * c.r - a2b1.i * c.i;
      coeffUFr[iDer][5].i = aux1.i + a2b1.r * c.i + a2b1.i * c.r;
      auxm1 = rho2 * (dpdb.r * uC2.r - dpdb.i * uC2.i);
      auxm2 = rho2 * (dpdb.r * uC2.i + dpdb.i * uC2.r);
      aux = coeffUFr[iDer][5].r * auxm1 - coeffUFr[iDer][5].i * auxm2;
      coeffUFr[iDer][5].i = coeffUFr[iDer][5].r * auxm2 + 
	                    coeffUFr[iDer][5].i * auxm1;
      /* Tsp */
      coeffUFr[iDer][5].r = aux; 
      
      /* computing the derivative */
      coeffUFr[iDer][5].r = (coeffUFr[iDer][5].r - coeffU[5].r) / DIV;
      coeffUFr[iDer][5].i = (coeffUFr[iDer][5].i - coeffU[5].i) / DIV;
      
      coeffUFr[iDer][6].r = aux1.r + a1b2.r * c.r - a1b2.i * c.i;
      coeffUFr[iDer][6].i = aux1.i + a1b2.r * c.i + a1b2.i * c.r;
      auxm1 = -rho2 * (dpda.r * uC2.r - dpda.i * uC2.i);
      auxm2 = -rho2 * (dpda.r * uC2.i + dpda.i * uC2.r);
      aux = coeffUFr[iDer][6].r * auxm1 - coeffUFr[iDer][6].i * auxm2;
      coeffUFr[iDer][6].i = coeffUFr[iDer][6].r * auxm2 + 
	                    coeffUFr[iDer][6].i * auxm1;          
      /* Tsp */
      coeffUFr[iDer][6].r = aux;
      
      /* computing the derivative */
      coeffUFr[iDer][6].r = (coeffUFr[iDer][6].r - coeffU[6].r) / DIV;
      coeffUFr[iDer][6].i = (coeffUFr[iDer][6].i - coeffU[6].i) / DIV;
      
      auxm1 = a1.r * aux2.r - a1.i * aux2.i;
      auxm2 = a1.r * aux2.i + a1.i * aux2.r;
      coeffUFr[iDer][7].r = auxm1 - (a2.r * aux3.r - a2.i * aux3.i);
      coeffUFr[iDer][7].i = auxm2 - (a2.r * aux3.i + a2.i * aux3.r);     
      auxm3 = 2 * rho2 * dpdb.r;
      auxm4 = 2 * rho2 * dpdb.i;
      aux = coeffUFr[iDer][7].r * auxm3 - coeffUFr[iDer][7].i * auxm4;
      coeffUFr[iDer][7].i = coeffUFr[iDer][7].r * auxm4 + 
	                    coeffUFr[iDer][7].i * auxm3;   
      /* Tss */
      coeffUFr[iDer][7].r = aux;
      
      /* computing the derivative */
      coeffUFr[iDer][7].r = (coeffUFr[iDer][7].r - coeffU[7].r) / DIV;
      coeffUFr[iDer][7].i = (coeffUFr[iDer][7].i - coeffU[7].i) / DIV;
   }
   /* restoring */
   vpFrechet = vpF;
   vsFrechet = vsF;
   rhoFrechet = rhoF;
}
/*                                                              */
/*  Function horSlownessFrechet()                               */
/*                                                              */
/*  Computing the following parameters used in the              */
/*  Frechet derivatives of the reflectivity matrices            */
/*                                                              */
/*  1) P-wave slowness squared                                  */
/*  2) S-wave slowness squared                                  */
/*  3) S-wave velocity squared                                  */
/*                                  ...for all layers           */
/*                                                              */
/*  Input parameters:                                           */
/*  alpha..................p-wave velocities of the model       */
/*                         global variable                      */
/*  beta...................s-wave velocities of the model       */
/*                         global variable                      */
/*  rho....................densities of the elastic model       */
/*                         global variable                      */
/*  thick..................thicknesses of the model             */
/*                         global variable                      */
/*  wC.....................complex frequency                    */
/*                         global variable                      */
/*  uC.....................complex slowness                     */
/*                         global variable                      */
/*                                                              */
/*  Output parameters:                                          */
/*  PSlowness[0..2][nL].......p-wave slowness squared           */
/*  SSlowness[0..2][nL].......s-wave slowness squared           */
/*  S2Velocity[0..2][nL]......s-wave velocity squared           */
/*                                                              */
/*  Notice that those quantities are also computed for the      */
/*  disturbed model, and are used in the derivatives of the     */
/*  reflection coefficients                                     */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
void horSlownessFrechet()
{
   /* declaration of variables */
   int iL, iF;                   /* counters */
   float cteQp1, cteQp2;
   float cteQs1, cteQs2;         /* constants used in absorption */
   float a, b;                   /* p and s-wave velocities */

   for (iF = 0; iF < 2; iF++)
   {
      for(iL = 0; iL <= nL; iL++)
      {
	 if (iF == 0)
	 {
	    a = alpha[iL];
	    b = beta[iL];
	 }
	 else
	 {
	    a = FACTOR * alpha[iL];
	    b = FACTOR * beta[iL];
	 }

	 /* constants used in absorption */
	 cteQp1 = a / (PI * qP[iL]);
	 cteQp2 = a / (2 * qP[iL]);
	 cteQs1 = b / (PI * qS[iL]);
	 cteQs2 = b / (2 * qS[iL]);

	 /* p-wave slowness squared */
	 PSlowness[iL][iF].r = a + cteQp1 * log(wCRwR);
	 PSlowness[iL][iF].i = cteQp2 - cteQp1 * wCP;

	 /* 1. / (PSlowness[iL] * PSlowness[iL]) */
	 auxm1 = PSlowness[iL][iF].r * PSlowness[iL][iF].r;
	 auxm2 = PSlowness[iL][iF].i * PSlowness[iL][iF].i;
	 aux = (auxm1 + auxm2) * (auxm1 + auxm2); aux = 1 / aux;
	 auxm3 = (auxm1 - auxm2) * aux;
	 PSlowness[iL][iF].i = -2 * PSlowness[iL][iF].r * 
	                            PSlowness[iL][iF].i * aux;
	 PSlowness[iL][iF].r = auxm3;
	 
	 /* computing an auxiliary quantity used in the gradient */
	 if (iF == 0)
	 {
	    derFactor[iL][0].r = PSlowness[iL][0].r / a;
	    derFactor[iL][0].i = PSlowness[iL][0].i / a;
	 }
      
	 /* s-wave velocity */
	 auxm1 = b + cteQs1 * log(wCRwR);
	 auxm2 = cteQs2 - cteQs1 * wCP;

	 /* S2Velocity[iL] * S2Velocity[iL] */
	 S2Velocity[iL][iF].r = auxm1 * auxm1 - auxm2 * auxm2;
	 S2Velocity[iL][iF].i = 2 * auxm1 * auxm2;

	 /* 1. / S2Velocity^2 */
	 aux = S2Velocity[iL][iF].r * S2Velocity[iL][iF].r + 
	       S2Velocity[iL][iF].i * S2Velocity[iL][iF].i;
	 SSlowness[iL][iF].r = S2Velocity[iL][iF].r / aux;
	 SSlowness[iL][iF].i = -S2Velocity[iL][iF].i / aux;

	 /* computing an auxiliary quantity used in the gradient */
	 if (iF == 0)
	 {
	    derFactor[iL][1].r = SSlowness[iL][0].r / b;
	    derFactor[iL][1].i = SSlowness[iL][0].i / b;
	 }

      }
   }
}