frechetslave.c

su 的源代码库

	       
		  /* irrI * (aux1, aux2) */
		  auxm1 = irrI[0][0].r * aux1.r - irrI[0][0].i * aux1.i;
		  auxm2 = irrI[0][0].r * aux1.i + irrI[0][0].i * aux1.r;
		  auxm3 = irrI[0][1].r * aux2.r - irrI[0][1].i * aux2.i;
		  auxm4 = irrI[0][1].r * aux2.i + irrI[0][1].i * aux2.r;
		  v2[iDer][0].r = auxm1 + auxm3;
		  v2[iDer][0].i = auxm2 + auxm4;
		  
		  auxm1 = irrI[1][0].r * aux1.r - irrI[1][0].i * aux1.i;
		  auxm2 = irrI[1][0].r * aux1.i + irrI[1][0].i * aux1.r;
		  auxm3 = irrI[1][1].r * aux2.r - irrI[1][1].i * aux2.i;
		  auxm4 = irrI[1][1].r * aux2.i + irrI[1][1].i * aux2.r;
		  v2[iDer][1].r = auxm1 + auxm3;
		  v2[iDer][1].i = auxm2 + auxm4;
	       }
	    }
	 
	    /* applying phase-shift to FRECHET derivatives */
	    /* loop over "active" layers */
	    for (iDer = 1; iDer <= numberPar * limRange; iDer++)
	    {
	       aux = v1[iDer][0].r * g[0].r - v1[iDer][0].i * g[0].i;
	       v1[iDer][0].i = v1[iDer][0].r * g[0].i + 
	 	               v1[iDer][0].i * g[0].r;
	       v1[iDer][0].r = aux;

	       aux = v1[iDer][1].r * g[1].r - v1[iDer][1].i * g[1].i;
	       v1[iDer][1].i = v1[iDer][1].r * g[1].i + 
	                       v1[iDer][1].i * g[1].r;
	       v1[iDer][1].r = aux;
	       
	       aux = v2[iDer][0].r * g[0].r - v2[iDer][0].i * g[0].i;
	       v2[iDer][0].i = v2[iDer][0].r * g[0].i + 
	                       v2[iDer][0].i * g[0].r;
	       v2[iDer][0].r = aux;
	    
	       aux = v2[iDer][1].r * g[1].r - v2[iDer][1].i * g[1].i;
	       v2[iDer][1].i = v2[iDer][1].r * g[1].i + 
		               v2[iDer][1].i * g[1].r;
	       v2[iDer][1].r = aux;
	    }
	    
	    /* compensating for free surface */
	    freeSurfaceFrechet(v1, v2);
	    
	    /* loop over offsets for computing the displacements */
	    displacementsFrechet(iU);
	 }
	 /* displacements in the radial or vertical direction */
	 /* (frequency domain) */
	 /* 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 "active" layers */
	 for (iDer = 0; iDer < numberPar * limRange; iDer++)
	 {
	    /* loop over offsets */
	    for (iR = 0; iR < info->nR; iR++)
	    {
	       /* radial ? */
	       if (RADIAL)
	       {
		  auxm1 = aux11[iDer][iR].r * dUCEp1.r - 
	 	          aux11[iDer][iR].i * dUCEp1.i;
		  auxm2 = aux11[iDer][iR].r * dUCEp1.i + 
		          aux11[iDer][iR].i * dUCEp1.r;
		  auxm3 = aux21[iDer][iR].r * dUCEp2.r - 
		          aux21[iDer][iR].i * dUCEp2.i;
		  auxm4 = aux21[iDer][iR].r * dUCEp2.i + 
		          aux21[iDer][iR].i * dUCEp2.r;
	 
		  displ.r = (auxm1 + auxm3) * wCCte.r - 
		            (auxm2 + auxm4) * wCCte.i;
		  displ.i = (auxm1 + auxm3) * wCCte.i + 
		            (auxm2 + auxm4) * wCCte.r;

		  /* filtering */
		  displ.r *= window[indexF] * SGN(recArray[iR]);
		  displ.i *= window[indexF] * SGN(recArray[iR]);
	       }
	       
	       if (VERTICAL)
	       {
		  auxm1 = aux12[iDer][iR].r * dUCEp1.r - 
		          aux12[iDer][iR].i * dUCEp1.i;
		  auxm2 = aux12[iDer][iR].r * dUCEp1.i + 
		          aux12[iDer][iR].i * dUCEp1.r;
		  auxm3 = aux22[iDer][iR].r * dUCEp2.r - 
		          aux22[iDer][iR].i * dUCEp2.i;
		  auxm4 = aux22[iDer][iR].r * dUCEp2.i + 
		          aux22[iDer][iR].i * dUCEp2.r;

		  displ.r = (auxm1 + auxm3) * wCCte.r - 
		            (auxm2 + auxm4) * wCCte.i;
		  displ.i = (auxm1 + auxm3) * wCCte.i + 
		            (auxm2 + auxm4) * wCCte.r;

		  /* filtering */
		  displ.r *= window[indexF];
		  displ.i *= window[indexF];
	       }
	    
	       /* computing the dot product of */
	       /* (residual^H . data covariance)  with the Frechet */
	       /* displacement vector */
	       grad[iDer] += resCD[iR][iF].r * displ.r + 
		             resCD[iR][iF].i * displ.i;
	       /* DD 
	       if (iDer == 0)
	       {
		  fprintf(logInfo.fp_log, "%f\n", grad[iDer]);
		  fprintf(logInfo.fp_log, "freq %f iR %d iDer %d grad %f\n", 
		  f, iR, iDer, grad[iDer]);
		  fprintf(logInfo.fp_log, 
			  "f %f res.r %f res.i %f displ.r %f displ.i %f\n", 
			  f,
			  resCD[iR][iF].r, resCD[iR][iF].i, displ.r, displ.i);
		  fflush(logInfo.fp_log);
	       }*/
	    }
	 }
      }
      /* sending partial gradient to the master process */
      SendFloat(grad, numberPar * limRange, masterId, PARTIAL_GRADIENT);
      
      /* DD 
      for (i = 0; i < numberPar * limRange; i++)
	 fprintf(logInfo.fp_log, "i %d grad %f\n", i, grad[i]);*/

      /* freeing memory */
      free2complex(resCD);

      /* 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();
}
/*                                                              */
/*  Function RmFrechet()                                        */
/*                                                              */
/*  Computing reflectivity matrices for a given model at a      */
/*  specific slowness and temporal frequency and its            */
/*  correspondent Frechet derivatives 	                        */
/*                                                              */
/*  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:                                          */
/*  rm.....................the reflectivity matrix              */
/*  DmB....................Frechet derivatives of the           */
/*                         reflectivity matrix                  */
/*  note the DmB data structure:                                */
/*  DmB[# of Layers][derivatives from the layer (or lim[0])     */
/*                   to the layer (or lim[1])][4 elements]      */
/*                                                              */
/*                                              Wences Gouveia  */
/*                                              September 1995  */
void RmFrechet()
{
   /* declaration of variables */
   int i, iL, iDer;             /* counter */
   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];

   /* checking depth limits */
   for (iDer = 0; iDer < numberPar; iDer++)
   {
      /* Vp -> Vs -> rho, if all active */
      if (lim[0] <= nL - 1 && nL - 1 < lim[1])
      {
	 /* Indexes: */
	 /*      [ 0  1 ]       */
	 /*      [ 2  3 ]       */
	 frechetRTd(nL - 1, nL, nL - 1);

	 DmB[nL - 1][limRange * iDer + 0][0] = coeffDFr[iDer][0];
	 DmB[nL - 1][limRange * iDer + 0][1] = coeffDFr[iDer][1];
	 DmB[nL - 1][limRange * iDer + 0][2] = coeffDFr[iDer][2];
	 DmB[nL - 1][limRange * iDer + 0][3] = coeffDFr[iDer][3];	       
      }
      if (lim[0] <= nL && nL < lim[1])
      {
	 frechetRTd(nL - 1, nL, nL);

 
	 DmB[nL - 1][limRange * iDer + 1][0] = coeffDFr[iDer][0];
	 DmB[nL - 1][limRange * iDer + 1][1] = coeffDFr[iDer][1];
	 DmB[nL - 1][limRange * iDer + 1][2] = coeffDFr[iDer][2];
	 DmB[nL - 1][limRange * iDer + 1][3] = coeffDFr[iDer][3];
      }
   }

   /* main loop over the nL layers */
   for (iL = nL - 1; iL >= 1; iL--)
   {
      /* square-root of PSlowness^2 - uuC */    
      auxm1 = PSlowness[iL][0].r - uuC.r;
      auxm2 = PSlowness[iL][0].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][0].r - uuC.r;
      auxm2 = SSlowness[iL][0].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);
      E[0][1].r = exp(auxm3) * cos(auxm4);
      E[0][1].i = exp(auxm3) * sin(auxm4);

      E[1][0] = E[0][1];

      /* cexp (bmI * wThick) */
      auxm1 = bmI.r * wThick.r - bmI.i * wThick.i;
      auxm2 = bmI.r * wThick.i + bmI.i * wThick.r;
      E[1][1].r = exp(auxm1) * cos(auxm2);
      E[1][1].i = exp(auxm1) * sin(auxm2);
	
      /* applying phase-shift */
      mT[0][0].r = mB[0][0].r * E[0][0].r - mB[0][0].i * E[0][0].i;
      mT[0][0].i = mB[0][0].r * E[0][0].i + mB[0][0].i * E[0][0].r;
      mT[0][1].r = mB[0][1].r * E[0][1].r - mB[0][1].i * E[0][1].i;
      mT[0][1].i = mB[0][1].r * E[0][1].i + mB[0][1].i * E[0][1].r;
      mT[1][0].r = mB[1][0].r * E[1][0].r - mB[1][0].i * E[1][0].i;
      mT[1][0].i = mB[1][0].r * E[1][0].i + mB[1][0].i * E[1][0].r;
      mT[1][1].r = mB[1][1].r * E[1][1].r - mB[1][1].i * E[1][1].i;
      mT[1][1].i = mB[1][1].r * E[1][1].i + mB[1][1].i * E[1][1].r;

      /* bottom-layer matrix - need a sequence of Ref and TRANS coeff. */
      RTd(iL - 1, iL);

      rD[0][0] = coeffD[0]; rD[0][1] = coeffD[1]; 
      rD[1][0] = coeffD[2]; rD[1][1] = coeffD[3];
      
      tD[0][0] = coeffD[4]; tD[0][1] = coeffD[5]; 
      tD[1][0] = coeffD[6]; tD[1][1] = coeffD[7];

      /* computing mT * tD */
      mTtD[0][0].r = mT[0][0].r * tD[0][0].r - mT[0][0].i * tD[0][0].i
	           + mT[0][1].r * tD[1][0].r - mT[0][1].i * tD[1][0].i;
      mTtD[0][0].i = mT[0][0].r * tD[0][0].i + mT[0][0].i * tD[0][0].r
	           + mT[0][1].r * tD[1][0].i + mT[0][1].i * tD[1][0].r;
      mTtD[0][1].r = mT[0][0].r * tD[0][1].r - mT[0][0].i * tD[0][1].i
	           + mT[0][1].r * tD[1][1].r - mT[0][1].i * tD[1][1].i;
      mTtD[0][1].i = mT[0][0].r * tD[0][1].i + mT[0][0].i * tD[0][1].r
	           + mT[0][1].r * tD[1][1].i + mT[0][1].i * tD[1][1].r;
      mTtD[1][0].r = mT[1][0].r * tD[0][0].r - mT[1][0].i * tD[0][0].i