modslave.c

su 的源代码库

   c.i = 2 * (rho[i1] * S2Velocity[i1].i - rho[i2] * S2Velocity[i2].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 - (rho[i1] - rho[i2]);
   aux1.i = cuu.i;

   aux2.r = cuu.r + rho[i2];
   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 - rho[i1];
   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;
   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) */
   coeffD[0].r = auxm1 * dd.r - auxm2 * dd.i;
   coeffD[0].i = auxm1 * dd.i + auxm2 * dd.r;                 /* Rpp */
   
   auxm1 = a2b2.r * caux3.r - a2b2.i * caux3.i;
   auxm2 = a2b2.r * caux3.i + a2b2.i * caux3.r;
   coeffD[1].r = aux1aux2.r + auxm1;
   coeffD[1].i = aux1aux2.i + auxm2;                          /* Rsp */  
								 
   coeffD[2].r = coeffD[1].r;
   coeffD[2].i = coeffD[1].i;                                 /* Rps*/  

   auxm3 = dpdb.r * uC2.r - dpdb.i * uC2.i;
   auxm4 = dpdb.r * uC2.i + dpdb.i * uC2.r;
   
   aux = auxm3 * coeffD[1].r - auxm4 * coeffD[1].i;
   coeffD[1].i = auxm3 * coeffD[1].i + auxm4 * coeffD[1].r;
   coeffD[1].r = aux;                                         /* Rsp */
  
   auxm3 = -dpda.r * uC2.r + dpda.i * uC2.i;
   auxm4 = -dpda.r * uC2.i - dpda.i * uC2.r;
   
   aux = auxm3 * coeffD[2].r - auxm4 * coeffD[2].i;
   coeffD[2].i = auxm3 * coeffD[2].i + auxm4 * coeffD[2].r;
   coeffD[2].r = aux;                                         /* Rps */

   auxm1 = d2d.r - d1d.r;
   auxm2 = d2d.i - d1d.i;
   auxm3 = 2 * rho1rho2 * (a1b2.r - a2b1.r);
   auxm4 = 2 * rho1rho2 * (a1b2.i - a2b1.i);

   coeffD[3].r = auxm1 - auxm3;
   coeffD[3].i = auxm2 - auxm4;
   aux = coeffD[3].r * dd.r - coeffD[3].i * dd.i;
   coeffD[3].i = coeffD[3].r * dd.i + coeffD[3].i * dd.r;
   coeffD[3].r = aux;                                         /* Rss */
  
   /* 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;
   coeffD[4].r = auxm1 - auxm3;     
   coeffD[4].i = auxm2 - auxm4;
   aux = 2 * rho[i1] * (dpda.r * coeffD[4].r - dpda.i * coeffD[4].i);
   coeffD[4].i = 2 * rho[i1] * (dpda.r * coeffD[4].i + dpda.i * coeffD[4].r);
   coeffD[4].r = aux;                                         /* Tpp */

   coeffD[5].r = aux1.r + a1b2.r * c.r - a1b2.i * c.i;
   coeffD[5].i = aux1.i + a1b2.r * c.i + a1b2.i * c.r;
   auxm1 = rho[i1] * (dpdb.r * uC2.r - dpdb.i * uC2.i);
   auxm2 = rho[i1] * (dpdb.r * uC2.i + dpdb.i * uC2.r);
   aux = coeffD[5].r * auxm1 - coeffD[5].i * auxm2;
   coeffD[5].i = coeffD[5].r * auxm2 + coeffD[5].i * auxm1;   /* Tsp */
   coeffD[5].r = aux; 

   coeffD[6].r = aux1.r + a2b1.r * c.r - a2b1.i * c.i;
   coeffD[6].i = aux1.i + a2b1.r * c.i + a2b1.i * c.r;
   auxm1 = -rho[i1] * (dpda.r * uC2.r - dpda.i * uC2.i);
   auxm2 = -rho[i1] * (dpda.r * uC2.i + dpda.i * uC2.r);
   aux = coeffD[6].r * auxm1 - coeffD[6].i * auxm2;
   coeffD[6].i = coeffD[6].r * auxm2 + coeffD[6].i * auxm1;   /* Tsp */
   coeffD[6].r = aux;

   auxm1 = a1.r * aux2.r - a1.i * aux2.i;
   auxm2 = a1.r * aux2.i + a1.i * aux2.r;
   coeffD[7].r = auxm1 - (a2.r * aux3.r - a2.i * aux3.i);
   coeffD[7].i = auxm2 - (a2.r * aux3.i + a2.i * aux3.r);     
   auxm3 = 2 * rho[i1] * dpdb.r;
   auxm4 = 2 * rho[i1] * dpdb.i;
   aux = coeffD[7].r * auxm3 - coeffD[7].i * auxm4;
   coeffD[7].i = coeffD[7].r * auxm4 + coeffD[7].i * auxm3;   /* Tss */
   coeffD[7].r = aux;
}

void RTu(int i1, int i2)
{
   /* declaration of variables */
   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 */
   
   /* square-root of PSlowness^2 - uuC */    
   auxm1 = PSlowness[i1].r - uuC.r;
   auxm2 = PSlowness[i1].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 = SSlowness[i1].r - uuC.r;
   auxm2 = SSlowness[i1].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 = PSlowness[i2].r - uuC.r;
   auxm2 = PSlowness[i2].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 = SSlowness[i2].r - uuC.r;
   auxm2 = SSlowness[i2].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);

   /* computing auxiliary quantities */
   rho1rho2 = rho[i1] * rho[i2];
   c.r = 2 * (rho[i1] * S2Velocity[i1].r - rho[i2] * S2Velocity[i2].r);
   c.i = 2 * (rho[i1] * S2Velocity[i1].i - rho[i2] * S2Velocity[i2].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 - (rho[i1] - rho[i2]);
   aux1.i = cuu.i;
   
   aux2.r = cuu.r + rho[i2];
   aux2.i = cuu.i;

   aux3.r = cuu.r - rho[i1];
   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) */
   coeffU[0].r = auxm1 * dd.r - auxm2 * dd.i;
   coeffU[0].i = auxm1 * dd.i + auxm2 * dd.r;                 /* Rpp */

   auxm1 = a1b1.r * caux2.r - a1b1.i * caux2.i;
   auxm2 = a1b1.r * caux2.i + a1b1.i * caux2.r;
   coeffU[1].r = aux1aux3.r + auxm1;
   coeffU[1].i = aux1aux3.i + auxm2;                          /* Rsp */  
								 
   coeffU[2].r = coeffU[1].r;
   coeffU[2].i = coeffU[1].i;                                 /* Rps*/  

   auxm3 = -dpdb.r * uC2.r + dpdb.i * uC2.i;
   auxm4 = -dpdb.r * uC2.i - dpdb.i * uC2.r;
   
   aux = auxm3 * coeffU[1].r - auxm4 * coeffU[1].i;
   coeffU[1].i = auxm3 * coeffU[1].i + auxm4 * coeffU[1].r;
   coeffU[1].r = aux;                                         /* Rsp */

   auxm3 = dpda.r * uC2.r - dpda.i * uC2.i;
   auxm4 = dpda.r * uC2.i + dpda.i * uC2.r;
   
   aux = auxm3 * coeffU[2].r - auxm4 * coeffU[2].i;
   coeffU[2].i = auxm3 * coeffU[2].i + auxm4 * coeffU[2].r;
   coeffU[2].r = aux;                                         /* Rps */

   auxm1 = d2u.r - d1u.r;
   auxm2 = d2u.i - d1u.i;
   auxm3 = 2 * rho1rho2 * (a2b1.r - a1b2.r);
   auxm4 = 2 * rho1rho2 * (a2b1.i - a1b2.i);

   coeffU[3].r = auxm1 - auxm3;
   coeffU[3].i = auxm2 - auxm4;
   aux = coeffU[3].r * dd.r - coeffU[3].i * dd.i;
   coeffU[3].i = coeffU[3].r * dd.i + coeffU[3].i * dd.r;
   coeffU[3].r = aux;                                         /* Rss */

   /* 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;
   coeffU[4].r = auxm1 - auxm3;     
   coeffU[4].i = auxm2 - auxm4;
   aux = 2 * rho[i2] * (dpda.r * coeffU[4].r - dpda.i * coeffU[4].i);
   coeffU[4].i = 2 * rho[i2] * (dpda.r * coeffU[4].i + dpda.i * coeffU[4].r);
   coeffU[4].r = aux;                                         /* Tpp */

   coeffU[5].r = aux1.r + a2b1.r * c.r - a2b1.i * c.i;
   coeffU[5].i = aux1.i + a2b1.r * c.i + a2b1.i * c.r;
   auxm1 = rho[i2] * (dpdb.r * uC2.r - dpdb.i * uC2.i);
   auxm2 = rho[i2] * (dpdb.r * uC2.i + dpdb.i * uC2.r);
   aux = coeffU[5].r * auxm1 - coeffU[5].i * auxm2;
   coeffU[5].i = coeffU[5].r * auxm2 + coeffU[5].i * auxm1;   /* Tsp */
   coeffU[5].r = aux; 

   coeffU[6].r = aux1.r + a1b2.r * c.r - a1b2.i * c.i;
   coeffU[6].i = aux1.i + a1b2.r * c.i + a1b2.i * c.r;
   auxm1 = -rho[i2] * (dpda.r * uC2.r - dpda.i * uC2.i);
   auxm2 = -rho[i2] * (dpda.r * uC2.i + dpda.i * uC2.r);
   aux = coeffU[6].r * auxm1 - coeffU[6].i * auxm2;
   coeffU[6].i = coeffU[6].r * auxm2 + coeffU[6].i * auxm1;   /* Tsp */
   coeffU[6].r = aux;

   auxm1 = a1.r * aux2.r - a1.i * aux2.i;
   auxm2 = a1.r * aux2.i + a1.i * aux2.r;
   coeffU[7].r = auxm1 - (a2.r * aux3.r - a2.i * aux3.i);
   coeffU[7].i = auxm2 - (a2.r * aux3.i + a2.i * aux3.r);     
   auxm3 = 2 * rho[i2] * dpdb.r;
   auxm4 = 2 * rho[i2] * dpdb.i;
   aux = coeffU[7].r * auxm3 - coeffU[7].i * auxm4;
   coeffU[7].i = coeffU[7].r * auxm4 + coeffU[7].i * auxm3;   /* Tss */
   coeffU[7].r = aux;
}

/* 
   Computation of bessel functions order 0 and 1.
   Extracted from Numerical Recipes.
*/

void Bessels(float x)
{
   register float xd3, x2, x4, x6, x8, x10, x12, sqrx;
   if (x <= 2.75) 
   {
      xd3 = x / 3;
      x2 = xd3 * xd3;
      x4 = x2 * x2;
      x6 = x2 * x4;
      x8 = x4 * x4;
      x10 = x4 * x6;
      x12 = x6 * x6;
      
      J00 = 1.-2.2499997*x2+1.2656208*x4
              -0.3163866*x6+0.0444479*x8
  	      -0.0039444*x10+0.00021*x12;
      
      J11 = (0.5-0.56249985*x2+0.21093573*x4
	        -0.03954289*x6+0.00443319*x8
	        -0.00031761*x10+0.00001109*x12)*x;
   }
   else
   {
      sqrx = sqrt(2. / (PI * x));
      J00 = sqrx * cos(x - 0.25 * PI);
      J11 = sqrx * cos(x - 0.75 * PI);
   }
}
/*
   Building windowing operator for convolution with seismograms
*/
void filter(float percW)
{
   /* declaration of variables */
   int i, iF, indexF;            /* counters */
   int wL;                       /* window length */

   /* memory allocation */
   window = alloc1float(nSamples / 2 + 1);
   /* reseting array */
   for (i = 0; i < nSamples / 2 + 1; i++)
      window[i] = 0;

   wL = percW * nF;   wL = 2 * wL - 1;
   /* building hanning window */
   for (i = 0, indexF = NINT(f1 / dF), iF = 0; iF < nF; iF++, indexF++)
   {
      if (hanningFlag)
	 window[indexF] =
	    .42 - .5 * cos(2 * PI * (float) iF / ((float) (nF - 1))) +
            .08 * cos(4 * PI * (float) iF / ((float) (nF - 1)));
      else
      {
	 window[indexF] = 1;
	 if (iF < (wL - 1) / 2)
	 {
	    window[indexF] = 
	       .42 - .5 * cos(2 * PI * (float) i / ((float) (wL - 1))) +
  	       .08 * cos(4 * PI * (float) i / ((float) (wL - 1)));
	    i++;
	 }
	 else if (iF > nF - (wL - 1) / 2)
	 {
	    i++;
	    window[indexF] =
	       .42 - .5 * cos(2 * PI * (float) i / ((float) (wL - 1))) +
  	       .08 * cos(4 * PI * (float) i / ((float) (wL - 1)));
	 }
      }
   }
}