grad.c

su 的源代码库

	    
	       auxm3 = rp[1][0].r * v1[0][0].r - rp[1][0].i * v1[0][0].i;
	       auxm4 = rp[1][0].r * v1[0][0].i + rp[1][0].i * v1[0][0].r;
	       auxm3 += rp[1][1].r * v1[0][1].r - rp[1][1].i * v1[0][1].i 
		     + Cs1.r;
	       auxm4 += rp[1][1].r * v1[0][1].i + rp[1][1].i * v1[0][1].r 
		     + Cs1.i;

	       /* DmB[0][active layers][0 1 2 3] * */
	       /*                      ((auxm1, auxm2), (auxm3, auxm4)) */
	       aux1.r = auxm1 * DmB[0][i * limRange + iL][0].r 
		      - auxm2 * DmB[0][i * limRange + iL][0].i 
 	 	      + auxm3 * DmB[0][i * limRange + iL][1].r 
		      - auxm4 * DmB[0][i * limRange + iL][1].i;
	       aux1.i = auxm1 * DmB[0][i * limRange + iL][0].i 
		      + auxm2 * DmB[0][i * limRange + iL][0].r
		      + auxm3 * DmB[0][i * limRange + iL][1].i 
		      + auxm4 * DmB[0][i * limRange + iL][1].r;
	       
	       aux2.r = auxm1 * DmB[0][i * limRange + iL][2].r 
		      - auxm2 * DmB[0][i * limRange + iL][2].i 
		      + auxm3 * DmB[0][i * limRange + iL][3].r 
		      - auxm4 * DmB[0][i * limRange + iL][3].i;
	       aux2.i = auxm1 * DmB[0][i * limRange + iL][2].i 
		      + auxm2 * DmB[0][i * limRange + iL][2].r  
		      + auxm3 * DmB[0][i * limRange + iL][3].i 
		      + auxm4 * DmB[0][i * limRange + iL][3].r;
	       
	       /* 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;
	       v1[iDer][0].r = auxm1 + auxm3;
	       v1[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;
	       v1[iDer][1].r = auxm1 + auxm3;
	       v1[iDer][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][0].r = auxm1 + auxm3;
	 v2[0][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[0][1].r = auxm1 + auxm3;
	 v2[0][1].i = auxm2 + auxm4;

	 /* loop over "active" layers */
	 for (iDer = 1, i = 0; i < numberPar; i++)
	 {
	    /* i = 0 -> Vp  */
	    /* i = 1 -> Vs  */
	    /* i = 2 -> rho */
	    for (iL = MIN(lim[0], 2); iL < MIN(lim[0], 2) + limRange; 
		 iL++, iDer++)
	    {
	       /* rp * [v2[0], v2[1]] + (As2, Bs2) */
	       auxm1 = rp[0][0].r * v2[0][0].r - rp[0][0].i * v2[0][0].i;
	       auxm2 = rp[0][0].r * v2[0][0].i + rp[0][0].i * v2[0][0].r;
	       auxm1 += rp[0][1].r * v2[0][1].r - rp[0][1].i * v2[0][1].i 
		     + As2.r;
	       auxm2 += rp[0][1].r * v2[0][1].i + rp[0][1].i * v2[0][1].r 
		     + As2.i;

	       auxm3 = rp[1][0].r * v2[0][0].r - rp[1][0].i * v2[0][0].i;
	       auxm4 = rp[1][0].r * v2[0][0].i + rp[1][0].i * v2[0][0].r;
	       auxm3 += rp[1][1].r * v2[0][1].r - rp[1][1].i * v2[0][1].i 
		     + Cs2.r;
	       auxm4 += rp[1][1].r * v2[0][1].i + rp[1][1].i * v2[0][1].r 
		     + Cs2.i;

	       /* DmB[0][active layers][0 1 2 3] * */
	       /*                      ((auxm1, auxm2), (auxm3, auxm4)) */
	       aux1.r = auxm1 * DmB[0][i * limRange + iL][0].r 
	 	      - auxm2 * DmB[0][i * limRange + iL][0].i 
		      + auxm3 * DmB[0][i * limRange + iL][1].r 
		      - auxm4 * DmB[0][i * limRange + iL][1].i;
	       aux1.i = auxm1 * DmB[0][i * limRange + iL][0].i 
		      + auxm2 * DmB[0][i * limRange + iL][0].r 
		      + auxm3 * DmB[0][i * limRange + iL][1].i 
		      + auxm4 * DmB[0][i * limRange + iL][1].r;
	       
	       aux2.r = auxm1 * DmB[0][i * limRange + iL][2].r 
		      - auxm2 * DmB[0][i * limRange + iL][2].i  
		      + auxm3 * DmB[0][i * limRange + iL][3].r 
		      - auxm4 * DmB[0][i * limRange + iL][3].i;
	       aux2.i = auxm1 * DmB[0][i * limRange + iL][2].i 
		      + auxm2 * DmB[0][i * limRange + iL][2].r 
		      + auxm3 * DmB[0][i * limRange + iL][3].i 
		      + auxm4 * DmB[0][i * limRange + iL][3].r;
	       
	       /* 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 < 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;
	 
	       dpl.i = (auxm1 + auxm3) * wCCte.r - (auxm2 + auxm4) * wCCte.i;
	       dpl.i = (auxm1 + auxm3) * wCCte.i + (auxm2 + auxm4) * wCCte.r;

	       /* filtering */
	       dpl.r *= window[indexF] * SGN(recArray[iR]);
	       dpl.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;

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

	       /* filtering */
	       dpl.r *= window[indexF];
	       dpl.i *= window[indexF];
	    }
	    
	    /* storing displacements in matrix displ */
	    displ[iDer][iR][indexF] = dpl;
	 }
      }
   }

   /* going to time domain and correctig for tau */
   for (iDer = 0; iDer < numberPar; iDer++)
   {
      for (iL = 0; iL < limRange; iL++)
      {
	 for (iR = 0; iR < nR; iR++)
 	 {
	    pfacr(1, nSamples, displ[iDer * limRange + iL][iR], buffer);
	   
	    /* correcting for tau */
	    for (iT = 0; iT < nSamples; iT++)
	    {
	       buffer[iT] *= exp(tau * iT * dt);
	    }
 
	    /* copying to operator F */
	    iT1 = NINT(t1 / dt);
	    for (iT = 0; iT < nDM; iT++)
	    {
	       if (IMPEDANCE && vpFrechet && iDer == 0)
	       {
		  F[iDer * limRange + iL][iR * nDM + iT] = 
  	          buffer[iT1 + iT] / rho[iL + lim[0]];
	       }
	       else if (IMPEDANCE && vsFrechet && (iDer == 0 || iDer == 1))
	       {
		  F[iDer * limRange + iL][iR * nDM + iT] = 
		  buffer[iT1 + iT] / rho[iL + lim[0]];
	       }
	       else if (IMPEDANCE && rhoFrechet && iDer == 2)
	       {
		  F[iDer * limRange + iL][iR * nDM + iT] =
		  - alpha[iL + lim[0]] * F[iL][iR * nDM + iT] 
		  - beta[iL + lim[0]] * F[iL + limRange][iR * nDM + iT] 
  	          + buffer[iT1 + iT];
	       }
	       else if (!IMPEDANCE)
	       {
		  F[iDer * limRange + iL][iR * nDM + iT] = buffer[iT1 + iT] ;
	       }
	    }
	 }
      }
   }
   
   /* if in the IMPEDANCE domain rearrange matrix F */
   if (IMPEDANCE)
   {
      if (rhoFrechet && !ipFrechet && !isFrechet)
      {
	 for (iL = 0; iL < limRange; iL++)
	 {
	    for (iR = 0; iR < nR; iR++)
	    {
	       for (iT = 0; iT < nDM; iT++)
	       {
		  F[iL][iR * nDM + iT] = F[iL + 2 * limRange][iR * nDM + iT];
	       }
	    }
	 }
      }
      else if (rhoFrechet && ipFrechet && !isFrechet)
      {
	 for (iL = 0; iL < limRange; iL++)
	 {
	    for (iR = 0; iR < nR; iR++)
	    {
	       for (iT = 0; iT < nDM; iT++)
	       {
		  F[iL + limRange][iR * nDM + iT] = 
		     F[iL + 2 * limRange][iR * nDM + iT];
	       }
	    }
	 }
      }
      else if (rhoFrechet && !ipFrechet && isFrechet)   
      {
	 for (iL = 0; iL < limRange; iL++)
	 {
	    for (iR = 0; iR < nR; iR++)
	    {
	       for (iT = 0; iT < nDM; iT++)
	       {
		  F[iL][iR * nDM + iT] = F[iL + limRange][iR * nDM + iT];
		  F[iL + limRange][iR * nDM + iT] = 
		     F[iL + 2 * limRange][iR * nDM + iT];
	       }
	    }
	 }
      }
   }
   
   /* freeing memory */
   free3complex(displ);
   free1float(buffer);
}