modslave.c

su 的源代码库

	    auxm3 = aux2.r * aux2.r - aux2.i * aux2.i;
	    auxm4 = 2 * aux2.r * aux2.i;
	    
	    aux3.r = auxm1 + auxm3;
	    aux3.i = auxm2 + auxm4;
	    
	    /* aux3 = 1. / aux3 */
	    aux = aux3.r * aux3.r + aux3.i * aux3.i;
	    aux3.r = aux3.r / aux;
	    aux3.i = -aux3.i / aux;
	    
	    /* ambm * aux1 * aux3 */
	    auxm1 = ambm.r * aux1.r - ambm.i * aux1.i;
	    auxm2 = ambm.r * aux1.i + ambm.i * aux1.r;
	    h[0][0].r = auxm1 * aux3.r - auxm2 * aux3.i;
	    h[0][0].i = auxm1 * aux3.i + auxm2 * aux3.r;
	    
	    /* aux2 * aux3 */
	    auxm1 = aux2.r * aux3.r - aux2.i * aux3.i;
	    auxm2 = aux2.r * aux3.i + aux2.i * aux3.r;
	    
	    /* bm * aux2 * aux3 */
	    h[0][1].r = auxm1 * bm.r - auxm2 * bm.i;
	    h[0][1].i = auxm1 * bm.i + auxm2 * bm.r;
	    
	    /* am * aux2 * aux3 */
	    h[1][0].r = auxm1 * am.r - auxm2 * am.i;
	    h[1][0].i = auxm1 * am.i + auxm2 * am.r;
	    
	    h[1][1].r = -h[0][0].r;
	    h[1][1].i = -h[0][0].i;
	    
	    /* computing phase shift (that's the matrix G in Muller's */
	    /* paper eq. (87) */
	    /* exp(j * am * wC * (-zs)) */
	    auxm1 = zs * (- amI.r * wC.r + amI.i * wC.i);
	    auxm2 = -zs * (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 * (-zs)) */
	    auxm1 = zs * (- bmI.r * wC.r + bmI.i * wC.i);
	    auxm2 = -zs * (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;
	    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;
	    }
	 }
	 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) */
	    if (RADIAL)
	    {
	       /* there's a 2 (free surface) / 2 (trapezoidal integration) */
	       /* simplified in the equation below */
	       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;
	    
	       uW[iR][iF].r = (auxm1 + auxm3) * wCCte.r - 
		  (auxm2 + auxm4) * wCCte.i;
	       uW[iR][iF].i = (auxm1 + auxm3) * wCCte.i + 
		  (auxm2 + auxm4) * wCCte.r;

	       /* filtering */
	       uW[iR][iF].r *= window[indexF] * SGN(recArray[iR]);
	       uW[iR][iF].i *= window[indexF] * SGN(recArray[iR]);
	    }

	    if (VERTICAL)
	    {
	       /* 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;
	       
	       uW[iR][iF].r = (auxm1 + auxm3) * wCCte.r - 
		  (auxm2 + auxm4) * wCCte.i;
	       uW[iR][iF].i = (auxm1 + auxm3) * wCCte.i + 
		  (auxm2 + auxm4) * wCCte.r;
	    
	       /* filtering */
	       uW[iR][iF].r *= window[indexF];
	       uW[iR][iF].i *= window[indexF];

	       /* DD 
	       fprintf(logInfo.fp_log, "iF %d fr %f fi %f window %f\n", 
		       iF, uW[iR][iF].r,uW[iR][iF].i, window[indexF]);
	       fflush(logInfo.fp_log);*/
	    }
	 }
      }
      /* sending frequency partition to the master process */
      SendCplx(uW[0], nFreqProc * nR, masterId, FREQUENCY_PARTITION);
      
      /* 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;
   }
}
/*
*/
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;