frechet.c

su 的源代码库

      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 = 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;
      d1d.r = auxm3 + a2b1.r * rho1rho2;
      d1d.i = auxm4 + a2b1.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;
      d1d.r += auxm3;
      d1d.i += auxm4;
      
      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;
      d2d.r = auxm3 + a1b2.r * rho1rho2;
      d2d.i = auxm4 + a1b2.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;
      d2d.r += auxm3;
      d2d.i += auxm4;
      
      /* more auxiliar quantities */
      /* d1d + d2d */
      dd.r = d1d.r + d2d.r;
      dd.i = d1d.i + d2d.i;
      
      /* 1 / dd */
      aux = dd.r * dd.r + dd.i * dd.i;
      dd.r = dd.r / aux;
      dd.i = -dd.i / aux;
      
      dpda.r = a1.r * dd.r - a1.i * dd.i;
      dpda.i = a1.r * dd.i + a1.i * dd.r;
      dpdb.r = b1.r * dd.r - b1.i * dd.i;
      dpdb.i = b1.r * dd.i + b1.i * dd.r;
      
      /* computing the coefficients - first reflection */
      auxm1 = d2d.r - d1d.r;
      auxm2 = d2d.i - d1d.i;
      /* (d2d - d1d) / (d1d + d2d) */
      coeffDFr[iDer][0].r = auxm1 * dd.r - auxm2 * dd.i;
      coeffDFr[iDer][0].i = auxm1 * dd.i + auxm2 * dd.r;     /* Rpp */

      /* computing the derivative */
      coeffDFr[iDer][0].r = (coeffDFr[iDer][0].r - coeffD[0].r) / DIV;
      coeffDFr[iDer][0].i = (coeffDFr[iDer][0].i - coeffD[0].i) / DIV;
      
      auxm1 = a2b2.r * caux3.r - a2b2.i * caux3.i;
      auxm2 = a2b2.r * caux3.i + a2b2.i * caux3.r;
      coeffDFr[iDer][1].r = aux1aux2.r + auxm1;
      coeffDFr[iDer][1].i = aux1aux2.i + auxm2;              /* Rsp */  
      
      coeffDFr[iDer][2].r = coeffDFr[iDer][1].r;
      coeffDFr[iDer][2].i = coeffDFr[iDer][1].i;             /* Rps */  
   
      auxm3 = dpdb.r * uC2.r - dpdb.i * uC2.i;
      auxm4 = dpdb.r * uC2.i + dpdb.i * uC2.r;
      
      aux = auxm3 * coeffDFr[iDer][1].r - auxm4 * coeffDFr[iDer][1].i;
      coeffDFr[iDer][1].i = auxm3 * coeffDFr[iDer][1].i + 
	                    auxm4 * coeffDFr[iDer][1].r;
      coeffDFr[iDer][1].r = aux;                             /* Rsp */
   
      /* computing the derivative */
      coeffDFr[iDer][1].r = (coeffDFr[iDer][1].r - coeffD[1].r) / DIV;
      coeffDFr[iDer][1].i = (coeffDFr[iDer][1].i - coeffD[1].i) / DIV;
      
      auxm3 = -dpda.r * uC2.r + dpda.i * uC2.i;
      auxm4 = -dpda.r * uC2.i - dpda.i * uC2.r;
      
      aux = auxm3 * coeffDFr[iDer][2].r - auxm4 * coeffDFr[iDer][2].i;
      coeffDFr[iDer][2].i = auxm3 * coeffDFr[iDer][2].i + 
	                    auxm4 * coeffDFr[iDer][2].r;
      coeffDFr[iDer][2].r = aux;                             /* Rps */
      
      /* computing the derivative */
      coeffDFr[iDer][2].r = (coeffDFr[iDer][2].r - coeffD[2].r) / DIV;
      coeffDFr[iDer][2].i = (coeffDFr[iDer][2].i - coeffD[2].i) / DIV;
      auxm1 = d2d.r - d1d.r;
      auxm2 = d2d.i - d1d.i;
      auxm3 = 2 * rho1rho2 * (a1b2.r - a2b1.r);
      auxm4 = 2 * rho1rho2 * (a1b2.i - a2b1.i);
      
      coeffDFr[iDer][3].r = auxm1 - auxm3;
      coeffDFr[iDer][3].i = auxm2 - auxm4;
      aux = coeffDFr[iDer][3].r * dd.r - coeffDFr[iDer][3].i * dd.i;
      coeffDFr[iDer][3].i = coeffDFr[iDer][3].r * dd.i + 
	                    coeffDFr[iDer][3].i * dd.r;
      coeffDFr[iDer][3].r = aux;                           /* Rss */
   
      /* computing the derivative */
      coeffDFr[iDer][3].r = (coeffDFr[iDer][3].r - coeffD[3].r) / DIV;
      coeffDFr[iDer][3].i = (coeffDFr[iDer][3].i - coeffD[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;
      coeffDFr[iDer][4].r = auxm1 - auxm3;     
      coeffDFr[iDer][4].i = auxm2 - auxm4;
      aux = 2 * rho1 * (dpda.r * coeffDFr[iDer][4].r - 
			dpda.i * coeffDFr[iDer][4].i);
      coeffDFr[iDer][4].i = 2 * rho1 * (dpda.r * coeffDFr[iDer][4].i + 
					dpda.i * coeffDFr[iDer][4].r);
      coeffDFr[iDer][4].r = aux;                             /* Tpp */
      
      /* computing the derivative */
      coeffDFr[iDer][4].r = (coeffDFr[iDer][4].r - coeffD[4].r) / DIV;
      coeffDFr[iDer][4].i = (coeffDFr[iDer][4].i - coeffD[4].i) / DIV;
      
      coeffDFr[iDer][5].r = aux1.r + a1b2.r * c.r - a1b2.i * c.i;
      coeffDFr[iDer][5].i = aux1.i + a1b2.r * c.i + a1b2.i * c.r;
      auxm1 = rho1 * (dpdb.r * uC2.r - dpdb.i * uC2.i);
      auxm2 = rho1 * (dpdb.r * uC2.i + dpdb.i * uC2.r);
      aux = coeffDFr[iDer][5].r * auxm1 - coeffDFr[iDer][5].i * auxm2;
      coeffDFr[iDer][5].i = coeffDFr[iDer][5].r * auxm2 + 
	                    coeffDFr[iDer][5].i * auxm1;   
      coeffDFr[iDer][5].r = aux;                             /* Tsp */
      
      /* computing the derivative */
      coeffDFr[iDer][5].r = (coeffDFr[iDer][5].r - coeffD[5].r) / DIV;
      coeffDFr[iDer][5].i = (coeffDFr[iDer][5].i - coeffD[5].i) / DIV;
      
      coeffDFr[iDer][6].r = aux1.r + a2b1.r * c.r - a2b1.i * c.i;
      coeffDFr[iDer][6].i = aux1.i + a2b1.r * c.i + a2b1.i * c.r;
      auxm1 = -rho1 * (dpda.r * uC2.r - dpda.i * uC2.i);
      auxm2 = -rho1 * (dpda.r * uC2.i + dpda.i * uC2.r);
      aux = coeffDFr[iDer][6].r * auxm1 - coeffDFr[iDer][6].i * auxm2;
      coeffDFr[iDer][6].i = coeffDFr[iDer][6].r * auxm2 + 
	                    coeffDFr[iDer][6].i * auxm1;  
      coeffDFr[iDer][6].r = aux;                             /* Tsp */
      
      /* computing the derivative */
      coeffDFr[iDer][6].r = (coeffDFr[iDer][6].r - coeffD[6].r) / DIV;
      coeffDFr[iDer][6].i = (coeffDFr[iDer][6].i - coeffD[6].i) / DIV;
      
      auxm1 = a1.r * aux2.r - a1.i * aux2.i;
      auxm2 = a1.r * aux2.i + a1.i * aux2.r;
      coeffDFr[iDer][7].r = auxm1 - (a2.r * aux3.r - a2.i * aux3.i);
      coeffDFr[iDer][7].i = auxm2 - (a2.r * aux3.i + a2.i * aux3.r);     
      auxm3 = 2 * rho1 * dpdb.r;
      auxm4 = 2 * rho1 * dpdb.i;
      aux = coeffDFr[iDer][7].r * auxm3 - coeffDFr[iDer][7].i * auxm4;
      coeffDFr[iDer][7].i = coeffDFr[iDer][7].r * auxm4 + 
	                    coeffDFr[iDer][7].i * auxm3; 
      coeffDFr[iDer][7].r = aux;                             /* Tss */
      
      /* computing the derivative */
      coeffDFr[iDer][7].r = (coeffDFr[iDer][7].r - coeffD[7].r) / DIV;
      coeffDFr[iDer][7].i = (coeffDFr[iDer][7].i - coeffD[7].i) / DIV;
   }
   /* restoring */
   vpFrechet = vpF;
   vsFrechet = vsF;
   rhoFrechet = rhoF;
}
/*                                                              */
/*  Function frechetRTu()                                       */
/*                                                              */
/*  Numerical computation of the derivatives of the             */
/*  reflection and transmission coefficients of upgoing         */
/*  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 te derivatives are with */
/*                         respect to the first or second layer */
/*                                                              */
/*  Output parameters:                                          */
/*  coeffUFr...............Derivatives of the elastic           */
/*                         reflection coefficients with respect */
/*                         to Vp, Vs and rho                    */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
void frechetRTu(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, caux2, aux1aux3;   /* auxiliar quantities */
   complex a1, a2, b1, b2;         /* all vertical slownesses */
   complex a1b1, a1b2, a2b1, a2b2; /* auxiliar quantities */
   complex a1a2b1b2;               /* auxiliar quantities */
   complex d1u, d2u;               /* 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;