sudrefsvslave.c

su 的源代码库

	    
	    /* computing phase shift (that's the matrix G in Muller's */
	    /* paper eq. (87) */
	    /* exp(j * am * wC * (-zm)) */
	    auxm1 = zm * (- amI.r * wC.r + amI.i * wC.i);
	    auxm2 = -zm * (amI.r * wC.i + amI.i * wC.r);
	    g[0].r = exp(auxm1) * cos(auxm2);
	    g[0].i = exp(auxm1) * sin(auxm2);
	    
	    /* exp(j * bm * wC * (-zm)) */
	    auxm1 = zm * (- bmI.r * wC.r + bmI.i * wC.i);
	    auxm2 = -zm * (bmI.r * wC.i + bmI.i * wC.r);
	    g[1].r = exp(auxm1) * cos(auxm2);
	    g[1].i = exp(auxm1) * sin(auxm2);
	    
	    /* computing the product I - R-R+ */
	    auxm1 = rm[0][0].r * rp[0][0].r - rm[0][0].i * rp[0][0].i;
	    auxm2 = rm[0][0].r * rp[0][0].i + rm[0][0].i * rp[0][0].r;
	    auxm3 = rm[0][1].r * rp[1][0].r - rm[0][1].i * rp[1][0].i;
	    auxm4 = rm[0][1].r * rp[1][0].i + rm[0][1].i * rp[1][0].r;
	    irr[0][0].r = 1 - (auxm1 + auxm3);
	    irr[0][0].i = - (auxm2 + auxm4);
	    
	    auxm1 = rm[0][0].r * rp[0][1].r - rm[0][0].i * rp[0][1].i;
	    auxm2 = rm[0][0].r * rp[0][1].i + rm[0][0].i * rp[0][1].r;
	    auxm3 = rm[0][1].r * rp[1][1].r - rm[0][1].i * rp[1][1].i;
	    auxm4 = rm[0][1].r * rp[1][1].i + rm[0][1].i * rp[1][1].r;
	    irr[0][1].r = - (auxm1 + auxm3);
	    irr[0][1].i = - (auxm2 + auxm4);
	    
	    auxm1 = rm[1][0].r * rp[0][0].r - rm[1][0].i * rp[0][0].i;
	    auxm2 = rm[1][0].r * rp[0][0].i + rm[1][0].i * rp[0][0].r;
	    auxm3 = rm[1][1].r * rp[1][0].r - rm[1][1].i * rp[1][0].i;
	    auxm4 = rm[1][1].r * rp[1][0].i + rm[1][1].i * rp[1][0].r;
	    irr[1][0].r = - (auxm1 + auxm3);
	    irr[1][0].i = - (auxm2 + auxm4);
	    
	    auxm1 = rm[1][0].r * rp[0][1].r - rm[1][0].i * rp[0][1].i;
	    auxm2 = rm[1][0].r * rp[0][1].i + rm[1][0].i * rp[0][1].r;
	    auxm3 = rm[1][1].r * rp[1][1].r - rm[1][1].i * rp[1][1].i;
	    auxm4 = rm[1][1].r * rp[1][1].i + rm[1][1].i * rp[1][1].r;
	    irr[1][1].r = 1 - (auxm1 + auxm3);
	    irr[1][1].i = - (auxm2 + auxm4);
	    
	    /* inverting irr explicitly */
	    auxm1 = irr[0][0].r * irr[1][1].r - irr[0][0].i * irr[1][1].i;
	    auxm2 = irr[0][0].r * irr[1][1].i + irr[0][0].i * irr[1][1].r;
	    auxm3 = irr[0][1].r * irr[1][0].r - irr[0][1].i * irr[1][0].i;
	    auxm4 = irr[0][1].r * irr[1][0].i + irr[0][1].i * irr[1][0].r;
	    aux1.r = auxm1 - auxm3;
	    aux1.i = auxm2 - auxm4;
	    
	    /* 1 / aux1 */
	    aux = aux1.r * aux1.r + aux1.i * aux1.i;
	    aux = 1. / aux;
	    aux1.r = aux1.r * aux;
	    aux1.i = -aux1.i * aux;
	    
	    /* Inverse of irr */
	    irrI[0][0].r = irr[1][1].r * aux1.r - irr[1][1].i * aux1.i;
	    irrI[0][0].i = irr[1][1].r * aux1.i + irr[1][1].i * aux1.r;
	    
	    irrI[0][1].r = -(irr[0][1].r * aux1.r - irr[0][1].i * aux1.i);
	    irrI[0][1].i = -(irr[0][1].r * aux1.i + irr[0][1].i * aux1.r);
	    
	    irrI[1][0].r = -(irr[1][0].r * aux1.r - irr[1][0].i * aux1.i);
	    irrI[1][0].i = -(irr[1][0].r * aux1.i + irr[1][0].i * aux1.r);
	    
	    irrI[1][1].r = irr[0][0].r * aux1.r - irr[0][0].i * aux1.i;
	    irrI[1][1].i = irr[0][0].r * aux1.i + irr[0][0].i * aux1.r;
	    
	    /* computing vectors V1,2, check eq (76) Muller's paper */
	    auxm1 = As1.r * rm[0][0].r - As1.i * rm[0][0].i;
	    auxm2 = As1.r * rm[0][0].i + As1.i * rm[0][0].r;
	    auxm3 = Cs1.r * rm[0][1].r - Cs1.i * rm[0][1].i;
	    auxm4 = Cs1.r * rm[0][1].i + Cs1.i * rm[0][1].r;
	    aux1.r = Bs1.r + (auxm1 + auxm3);
	    aux1.i = Bs1.i + (auxm2 + auxm4);
	    
	    auxm1 = As1.r * rm[1][0].r - As1.i * rm[1][0].i;
	    auxm2 = As1.r * rm[1][0].i + As1.i * rm[1][0].r;
	    auxm3 = Cs1.r * rm[1][1].r - Cs1.i * rm[1][1].i;
	    auxm4 = Cs1.r * rm[1][1].i + Cs1.i * rm[1][1].r;
	    aux2.r = Ds1.r + (auxm1 + auxm3);
	    aux2.i = Ds1.i + (auxm2 + auxm4);
	    
	    auxm1 = aux1.r * irrI[0][0].r - aux1.i * irrI[0][0].i;
	    auxm2 = aux1.r * irrI[0][0].i + aux1.i * irrI[0][0].r;
	    auxm3 = aux2.r * irrI[0][1].r - aux2.i * irrI[0][1].i;
	    auxm4 = aux2.r * irrI[0][1].i + aux2.i * irrI[0][1].r;
	    v1[0].r = auxm1 + auxm3;
	    v1[0].i = auxm2 + auxm4;
	    
	    auxm1 = aux1.r * irrI[1][0].r - aux1.i * irrI[1][0].i;
	    auxm2 = aux1.r * irrI[1][0].i + aux1.i * irrI[1][0].r;
	    auxm3 = aux2.r * irrI[1][1].r - aux2.i * irrI[1][1].i;
	    auxm4 = aux2.r * irrI[1][1].i + aux2.i * irrI[1][1].r;
	    v1[1].r = auxm1 + auxm3;
	    v1[1].i = auxm2 + auxm4;
	    
	    auxm1 = As2.r * rm[0][0].r - As2.i * rm[0][0].i;
	    auxm2 = As2.r * rm[0][0].i + As2.i * rm[0][0].r;
	    auxm3 = Cs2.r * rm[0][1].r - Cs2.i * rm[0][1].i;
	    auxm4 = Cs2.r * rm[0][1].i + Cs2.i * rm[0][1].r;
	    aux1.r = Bs2.r + (auxm1 + auxm3);
	    aux1.i = Bs2.i + (auxm2 + auxm4);

	    auxm1 = As2.r * rm[1][0].r - As2.i * rm[1][0].i;
	    auxm2 = As2.r * rm[1][0].i + As2.i * rm[1][0].r;
	    auxm3 = Cs2.r * rm[1][1].r - Cs2.i * rm[1][1].i;
	    auxm4 = Cs2.r * rm[1][1].i + Cs2.i * rm[1][1].r;
	    aux2.r = Ds2.r + (auxm1 + auxm3);
	    aux2.i = Ds2.i + (auxm2 + auxm4);
	    
	    auxm1 = aux1.r * irrI[0][0].r - aux1.i * irrI[0][0].i;
	    auxm2 = aux1.r * irrI[0][0].i + aux1.i * irrI[0][0].r;
	    auxm3 = aux2.r * irrI[0][1].r - aux2.i * irrI[0][1].i;
	    auxm4 = aux2.r * irrI[0][1].i + aux2.i * irrI[0][1].r;
	    v2[0].r = auxm1 + auxm3;
	    v2[0].i = auxm2 + auxm4;
	    
	    auxm1 = aux1.r * irrI[1][0].r - aux1.i * irrI[1][0].i;
	    auxm2 = aux1.r * irrI[1][0].i + aux1.i * irrI[1][0].r;
	    auxm3 = aux2.r * irrI[1][1].r - aux2.i * irrI[1][1].i;
	    auxm4 = aux2.r * irrI[1][1].i + aux2.i * irrI[1][1].r;
	    v2[1].r = auxm1 + auxm3;
	    v2[1].i = auxm2 + auxm4;

	    /* applying phase-shift */
	    aux = v1[0].r * g[0].r - v1[0].i * g[0].i;
	    v1[0].i = v1[0].r * g[0].i + v1[0].i * g[0].r;
	    v1[0].r = aux;
	    
	    aux = v1[1].r * g[1].r - v1[1].i * g[1].i;
	    v1[1].i = v1[1].r * g[1].i + v1[1].i * g[1].r;
	    v1[1].r = aux;
	    
	    aux = v2[0].r * g[0].r - v2[0].i * g[0].i;
	    v2[0].i = v2[0].r * g[0].i + v2[0].i * g[0].r;
	    v2[0].r = aux;
	    
	    aux = v2[1].r * g[1].r - v2[1].i * g[1].i;
	    v2[1].i = v2[1].r * g[1].i + v2[1].i * g[1].r;
	    v2[1].r = aux;
	    
	    /* multiplication by matrix h */
	    auxm1 = h[0][0].r * v1[0].r - h[0][0].i * v1[0].i;
	    auxm2 = h[0][0].r * v1[0].i + h[0][0].i * v1[0].r;
	    auxm3 = h[0][1].r * v1[1].r - h[0][1].i * v1[1].i;
	    auxm4 = h[0][1].r * v1[1].i + h[0][1].i * v1[1].r;
	    v1A = v1[0];              /* for future use */
	    v1[0].r = auxm1 + auxm3;
	    v1[0].i = auxm2 + auxm4;
	    
	    auxm1 = h[0][0].r * v2[0].r - h[0][0].i * v2[0].i;
	    auxm2 = h[0][0].r * v2[0].i + h[0][0].i * v2[0].r;
	    auxm3 = h[0][1].r * v2[1].r - h[0][1].i * v2[1].i;
	    auxm4 = h[0][1].r * v2[1].i + h[0][1].i * v2[1].r;
	    v2A = v2[0];              /* for future use */
	    v2[0].r = auxm1 + auxm3;
	    v2[0].i = auxm2 + auxm4;
	    
	    auxm1 = h[1][0].r * v1A.r - h[1][0].i * v1A.i;
	    auxm2 = h[1][0].r * v1A.i + h[1][0].i * v1A.r;
	    auxm3 = h[1][1].r * v1[1].r - h[1][1].i * v1[1].i;
	    auxm4 = h[1][1].r * v1[1].i + h[1][1].i * v1[1].r;
	    v1[1].r = auxm1 + auxm3;
	    v1[1].i = auxm2 + auxm4;
	    
	    auxm1 = h[1][0].r * v2A.r - h[1][0].i * v2A.i;
	    auxm2 = h[1][0].r * v2A.i + h[1][0].i * v2A.r;
	    auxm3 = h[1][1].r * v2[1].r - h[1][1].i * v2[1].i;
	    auxm4 = h[1][1].r * v2[1].i + h[1][1].i * v2[1].r;
	    v2[1].r = auxm1 + auxm3;
	    v2[1].i = auxm2 + auxm4;
	    
	    /* loop over offsets for computing the displacements */
	    for (iR = 0; iR < nR; iR++)
	    {
	       /* Bessel functions */
	       arg = sqrt((uC.r * wC.r - uC.i * wC.i) * 
			  (uC.r * wC.r - uC.i * wC.i) + 
			  (wC.i * uC.r + uC.i * wC.r) * 
			  (wC.i * uC.r + uC.i * wC.r)) * ABS(recArray[iR]);
	       Bessels(arg);
	       
	       /* r component */
	       aux1.r = -J11 * taper[iU] * v1[0].r;
	       aux1.i = -J11 * taper[iU] * v1[0].i;
	       
	       aux11[iR].r += aux1.r + aux11Old[iR].r;
	       aux11[iR].i += aux1.i + aux11Old[iR].i;
	       aux11Old[iR] = aux1;
	       
	       aux1.r = J00 * taper[iU] * v2[0].r;
	       aux1.i = J00 * taper[iU] * v2[0].i;
	       
	       aux21[iR].r += aux1.r + aux21Old[iR].r;
	       aux21[iR].i += aux1.i + aux21Old[iR].i;
	       aux21Old[iR] = aux1;
	       
	       /* z component */
	       aux1.r = -v1[1].i;
	       aux1.i = v1[1].r;
	       
	       aux1.r = J00 * taper[iU] * aux1.r;
	       aux1.i = J00 * taper[iU] * aux1.i;
	       
	       aux12[iR].r += aux1.r + aux12Old[iR].r;
	       aux12[iR].i += aux1.i + aux12Old[iR].i;
	       aux12Old[iR] = aux1;	
	       
	       aux1.r = -v2[1].i;
	       aux1.i =  v2[1].r;
	       
	       aux1.r = J11 * taper[iU] * aux1.r;
	       aux1.i = J11 * taper[iU] * aux1.i;
	       
	       aux22[iR].r += aux1.r + aux22Old[iR].r;
	       aux22[iR].i += aux1.i + aux22Old[iR].i;
	       aux22Old[iR] = aux1;
	    }
	 }
	 /* there's a 2 (free surface) / 2 (trapezoidal integration) */
	 /* simplified in the equation below */
	 dUCEp1.r = epslon1 * dUC.r;
	 dUCEp1.i = epslon1 * dUC.i;
	 dUCEp2.r = epslon2 * dUC.r;
	 dUCEp2.i = epslon2 * dUC.i;
	 
	 /* loop over offsets */
	 for (iR = 0; iR < nR; iR++)
	 {
	    
	    /* displacements in the radial direction (frequency domain) */
	    auxm1 = aux11[iR].r * dUCEp1.r - aux11[iR].i * dUCEp1.i;
	    auxm2 = aux11[iR].r * dUCEp1.i + aux11[iR].i * dUCEp1.r;
	    auxm3 = aux21[iR].r * dUCEp2.r - aux21[iR].i * dUCEp2.i;
	    auxm4 = aux21[iR].r * dUCEp2.i + aux21[iR].i * dUCEp2.r;
	    
	    uRw[iR][iF].r = (auxm1 + auxm3) * wCCte.r - 
	       (auxm2 + auxm4) * wCCte.i;
	    uRw[iR][iF].i = (auxm1 + auxm3) * wCCte.i + 
	       (auxm2 + auxm4) * wCCte.r;

	    /* filtering */
	    uRw[iR][iF].r *= window[indexF] * SGN(recArray[iR]);
	    uRw[iR][iF].i *= window[indexF] * SGN(recArray[iR]);
	    
	    /* displacements in the z direction (frequency domain) */
	    /* there's a 2 (free surface) / 2 (trapezoidal integration) */
	    /* simplified in the equation below */
	    auxm1 = aux12[iR].r * dUCEp1.r - aux12[iR].i * dUCEp1.i;
	    auxm2 = aux12[iR].r * dUCEp1.i + aux12[iR].i * dUCEp1.r;
	    auxm3 = aux22[iR].r * dUCEp2.r - aux22[iR].i * dUCEp2.i;
	    auxm4 = aux22[iR].r * dUCEp2.i + aux22[iR].i * dUCEp2.r;
	    
	    uZw[iR][iF].r = (auxm1 + auxm3) * wCCte.r - 
	       (auxm2 + auxm4) * wCCte.i;
	    uZw[iR][iF].i = (auxm1 + auxm3) * wCCte.i + 
	       (auxm2 + auxm4) * wCCte.r;

	    /* filtering */
	    uZw[iR][iF].r *= window[indexF];
	    uZw[iR][iF].i *= window[indexF];
	 }
      }
      /* sending frequency partition to the master process */
      SendCplx(uZw[0], nFreqProc * nR, masterId, FREQUENCY_PARTITION_VERTICAL);
      SendCplx(uRw[0], nFreqProc * nR, masterId, FREQUENCY_PARTITION_RADIAL);
      
      /* keep working ? */
      masterId = RecvInt(&die, 1, -1, DIE);
      
      if (verbose)
      {
	 fprintf(logInfo.fp_log, "Slave receiving flag to stop : %d\n", die);
	 fflush(logInfo.fp_log);
      }
      if (die) break;
   } while (FOREVER);

   /* quitting PVM */
   EndOfSlave();    
}

/*
   Computing frequency-dependent slowness (squared) for all layer of the
   model
*/
void horSlowness()
{
   /* declaration of variables */
   int iL;                       /* counter */
   float cteQp1, cteQp2;
   float cteQs1, cteQs2;         /* constants used in absorption */

   for(iL = 0; iL <= nL; iL++)
   {

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

      /* p-wave slowness squared */
      PSlowness[iL].r = alpha[iL] + cteQp1 * log(wCRwR);
      PSlowness[iL].i = cteQp2 - cteQp1 * wCP;
      
      /* 1. / (PSlowness[iL] * PSlowness[iL]) */
      auxm1 = PSlowness[iL].r * PSlowness[iL].r;
      auxm2 = PSlowness[iL].i * PSlowness[iL].i;
      aux = (auxm1 + auxm2) * (auxm1 + auxm2); aux = 1 / aux;
      auxm3 = (auxm1 - auxm2) * aux;
      PSlowness[iL].i = -2 * PSlowness[iL].r * PSlowness[iL].i * aux;
      PSlowness[iL].r = auxm3;

      /* s-wave velocity */
      auxm1 = beta[iL] + cteQs1 * log(wCRwR);
      auxm2 = cteQs2 - cteQs1 * wCP;

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

      /* 1. / S2Velocity^2 */
      aux = S2Velocity[iL].r * S2Velocity[iL].r + 
	    S2Velocity[iL].i * S2Velocity[iL].i;
      SSlowness[iL].r = S2Velocity[iL].r / aux;
      SSlowness[iL].i = -S2Velocity[iL].i / aux;
   }
}
/*
   The reflectivity matrices according to the formultation of G. Muller 
   (Journal of Geophysics, 1985, vol. 58) are implemented in this program.
   The elastic case is handled, and this routine makes calls to subroutines
   that compute the reflection and transmission coefficients for up and
   down going waves.


   Input parameters:
   




   Output parameters:




                                                   Wences - CWP
						   06-08-95
*/
void Rm()
{
   /* declaration of variables */
   int iL;                      /* layer index */
                                /* used in absorption */
   complex aux1;                /* auxiliar quantities */
   complex am, amI, bm, bmI;    /* vertical slownesses for P and S waves */
                                /* amI = am * I, bmI = bm * I */
   complex wThick;              /* wC * thickness */
   complex E[2][2];             /* phase-shift matrix */
   complex mT[2][2];            /* reflectivity matrix at top of layer */
   complex mTtD[2][2];          /* mt * tD */
   complex tUinv[2][2];         /* tU * inv */
   complex mB[2][2];            /* reflectivity matrix at bottom of layer */
   complex rD[2][2], tD[2][2];  /* reflec. and transm. coefficients for */
                                /* downgoing waves */
   complex rU[2][2], tU[2][2];  /* reflec. and transm. coefficients for */
                                /* upgoing waves */
   complex mAux[2][2];          /* auxiliar matrix */
   complex inv[2][2];           /* inv = (I - mT * rU)^-1 */
   
   /* initializing the reflectivity matrix at the bottom of half space */
   mT[0][0] = zeroC; mT[0][1] = zeroC;
   mT[1][0] = zeroC; mT[1][1] = zeroC;
   
   /* initializing the reflectivity matrix at the bottom */
   /* of layer nL - 1 */
   RTd(nL - 1, nL);

   mB[0][0] = coeffD[0]; mB[0][1] = coeffD[1]; 
   mB[1][0] = coeffD[2]; mB[1][1] = coeffD[3];

   /* main loop over the nL layers */
   for (iL = nL - 1; iL >= 1; iL--)
   {
      /* square-root of PSlowness^2 - uuC */    
      auxm1 = PSlowness[iL].r - uuC.r;
      auxm2 = PSlowness[iL].i - uuC.i;
      auxm3 = sqrt(auxm1 * auxm1 + auxm2 * auxm2);
      auxm3 = sqrt(auxm3);
      angle = atan2(auxm2, auxm1) / 2;
      am.r = auxm3 * cos(angle);
      am.i = auxm3 * sin(angle);

      /* am * I */
      amI.r = -am.i;
      amI.i = am.r;

      /* square-root of SSlowness^2 - uuC */    
      auxm1 = SSlowness[iL].r - uuC.r;
      auxm2 = SSlowness[iL].i - uuC.i;
      auxm3 = sqrt(auxm1 * auxm1 + auxm2 * auxm2);
      auxm3 = sqrt(auxm3);
      angle = atan2(auxm2, auxm1) / 2;
      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[iL]);
      wThick.i = wC.i * (-2 * thick[iL]);

      /* 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);