modslave.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

#include "modSlave.h"

void
main ()
{
   /* declaration of variables */
   
   FILE *fp;                     /* file pointer for receiver array */
   char *recFile = " ";          /* file with receiver coordinates */
   int indexF, iU, iF, iR, i;    /* counters */
   int masterId;                 /* master id */
   int die;                      /* if die=1 process stops */
   int nU;                       /* number of slownesses */
   int nR;                       /* number of receivers */
   int wL;                       /* taper length */
   int wInfo[2];                 /* frequency limits */
   int nFreqProc;                /* number of frequencies per processor */
   int directWave;               /* direct wave flag */
   int VERTICAL, RADIAL;         /* specifies geophone orientation */
   int verbose;                  /* dialogue flag */
   float tau;                    /* magnitude of wrap-around attenuation */
   float percU;                  /* percentual of slowness */
                                 /* domain that are windowed */
   float percW;                  /* percentual of frequency */
                                 /* domain that are windowed */
   float r1, dR;                 /* defines receiver array */
   float u1, u2;                 /* slowness window */
   float u, w, f;                /* slowness & frequency */
   float phi;                    /* azimuth angle */
   float F1, F2, F3;             /* source components */
   float epslon1, epslon2;       /* auxiliary quantities */ 
   float cte;                    /* = 4 PI rho */
   float *taper;                 /* taper for slowness domain */
   float arg;                    /* argument of the Bessel function */
   float dU;                     /* slowness interval */
   float wRef;                   /* used in complex slowness */
   float *recArray;              /* array with receiver coordinates */
   complex dUC;                  /* complex slowness interval */
   complex v1A, v2A;             /* auxiliary values */
   complex dUCEp1, dUCEp2;       /* dUC * epslon1 and dUC * epslon2 */
   complex muC;                  /* uC * -1 */
   complex wCCte;                /* auxiliar variable */
   complex am;                   /* vertical P-wave slownesses */
   complex amInv;                /* 1. / am */
   complex amI;                  /* amI = am * I */
   complex bm;                   /* vertical S-wave slownesses */
   complex bmInv;                /* 1. / bm */
   complex bmI;                  /* bmI = bm * I */
   complex ambm;                 /* - am * bm */
   complex As1, As2;             /* amplitudes of plane wave components (P)*/
   complex Cs1, Cs2;             /* amplitudes of plane wave components (S)*/
                                 /* downgoing waves */
   complex Bs1, Bs2;             /* amplitudes of plane wave components (P)*/
   complex Ds1, Ds2;             /* amplitudes of plane wave components (S)*/
                                 /* upgoing waves */
   complex irr[2][2];            /* = (I - R-R+) */
   complex irrI[2][2];           /* = (I - R-R+)^-1 */
   complex v1[2], v2[2];         /* potential vectors */
   complex g[2];                 /* phase-shift vector */
   complex h[2][2];              /* free-surface matrix */
   complex **uW;                 /* displacements in the vertical or radial */
                                 /* direction in the frequency domain */
   complex *aux11, *aux12, *aux21, *aux22;
   complex *aux11Old, *aux12Old, *aux21Old, *aux22Old;
   register complex aux1, aux2, aux3;     /* auxiliar quantities */
   INFO info[1];                 /* basic information for slaves */
   SeisSlave logInfo;            /* report structure */     

   /* logging information */
   InitLog(&logInfo, PROCESS_MODELING);

   /* receiving INFO data structure */
   masterId = RecvINFO(info, 1, -1, GENERAL_INFORMATION);

   recFile = info->recFile;
   directWave = info->directWave;
   nL = info->nL;
   nF = info->nF;
   nR = info->nR;
   nU = info->nU;
   hanningFlag = info->hanningFlag;
   nSamples = info->nSamples;
   nFreqProc = info->nFreqProc;
   VERTICAL = info->vertical;
   RADIAL = info->radial;
   r1 = info->r1;
   dR = info->dR;
   zs = info->zs;
   u1 = info->u1;
   u2 = info->u2;
   dU = info->dU;
   f1 = info->f1;
   f2 = info->f2;
   dF = info->dF;
   F1 = info->F1;
   F2 = info->F2;
   F3 = info->F3;
   percU = info->percU;
   percW = info->percW;
   wR = info->wR;
   tau = info->tau;
   
   if (verbose)
   {
      fprintf(logInfo.fp_log, "Starting modeling for %s\n", info->id);
      fflush(logInfo.fp_log);
   }

   /* constants */
   zeroC = cmplx(0,0);

   /* memory allocation */
   thick = alloc1float(nL + 1);
   alpha = alloc1float(nL + 1);
   beta = alloc1float(nL + 1);
   rho = alloc1float(nL + 1);
   qP = alloc1float(nL + 1);
   qS = alloc1float(nL + 1);
   PSlowness = alloc1complex(nL + 1);
   SSlowness = alloc1complex(nL + 1);
   S2Velocity = alloc1complex(nL + 1);
   recArray = alloc1float(nR);
   aux11 = alloc1complex(nR);
   aux12 = alloc1complex(nR);
   aux21 = alloc1complex(nR);
   aux22 = alloc1complex(nR);
   aux11Old = alloc1complex(nR);
   aux12Old = alloc1complex(nR);
   aux21Old = alloc1complex(nR);
   aux22Old = alloc1complex(nR);

   uW = alloc2complex(nFreqProc, nR);

   /* receiving elastic model */
   RecvFloat(thick, nL + 1, -1, THICKNESS);
   RecvFloat(rho, nL + 1, -1, DENSITY);
   RecvFloat(alpha, nL + 1, -1, ALPHAS);
   RecvFloat(qP, nL + 1, -1, QALPHA);
   RecvFloat(beta, nL + 1, -1, BETAS);
   RecvFloat(qS, nL + 1, -1, QBETA);

   /* DD 
   fprintf(logInfo.fp_log, "nF %d\n", nF);
   fprintf(logInfo.fp_log, "nFreqProc %d\n", nFreqProc);
   fprintf(logInfo.fp_log, "tau %f\n", tau);
   fprintf(logInfo.fp_log, "percW %f\n", percW);
   fprintf(logInfo.fp_log, "percU %f\n", percU);
   fprintf(logInfo.fp_log, "u1 %f\n", u1);
   fprintf(logInfo.fp_log, "u2 %f\n", u2);
   fprintf(logInfo.fp_log, "r1 %f\n", r1);
   fprintf(logInfo.fp_log, "dR %f\n", dR);
   fprintf(logInfo.fp_log, "F1 %f\n", F1);
   fprintf(logInfo.fp_log, "F2 %f\n", F2);
   fprintf(logInfo.fp_log, "F3 %f\n", F3);
   fprintf(logInfo.fp_log, "zs %f\n", zs);
   fprintf(logInfo.fp_log, "nSamples %d\n", nSamples);
   for (iU = 0; iU < nL; iU++)
   {
      fprintf(logInfo.fp_log, "t %f d %f a %f qa %f b %f qb %f\n", 
	      thick[iU], rho[iU], alpha[iU], qP[iU], beta[iU], qS[iU]);
   }
   fflush(logInfo.fp_log);*/

   /* source is at 1st layer */
   thick[0] -= zs;
   
   /* computing the window length for the slowness domain */
   epslon1 = (u2 - u1) * percU;
   wL = NINT(epslon1 / dU);
   wL = 2 * wL + 1;
   u2 += epslon1;
   nU = NINT((u2 - u1) / dU);    /* new nU to preserve last slowness */
                                 /* w/o being windowed */
   /* DD 
   fprintf(logInfo.fp_log, "u1 %f\n", u1);
   fprintf(logInfo.fp_log, "u2 %f\n", u2);
   fprintf(logInfo.fp_log, "dU %f nU %d\n", dU, nU);
   fflush(logInfo.fp_log);*/
      
   /* building window for slowness integration */
   taper = alloc1float(nU);
   for (i = (wL - 1) / 2, iU = 0; iU < nU; iU++)
   {
      taper[iU] = 1;
      if (iU >= nU - (wL - 1) / 2)
      {
         i++;
	 taper[iU] =
	    .42 - .5 * cos(2 * PI * (float) i / ((float) (wL - 1))) +
            .08 * cos(4 * PI * (float) i / ((float) (wL - 1)));
      }
   }

   /* building windowing operator */
   /* DD
   fprintf(logInfo.fp_log, "going to filter\n");
   fflush(logInfo.fp_log);*/
   filter(percW);
   
   /* I will assume that the receivers are in line (at z = 0) so phi = 0 */
   phi = 0;
   epslon1 = F3;
   epslon2 = F1 * cos(phi) + F2 * sin(phi);

   /* auxiliar constant */
   cte = 1. / (4 * PI * rho[0]);

   /* normalization for the complex slowness */
   if (f1 > 7.5)
      wRef = f1 * 2 * PI;
   else
      wRef = 7.5 * 2 * PI;

   /* specific geometry */
   fp = fopen(recFile, "r");
   if (fp == NULL)
   {
      /* standard end-on */
      for (iR = 0; iR < nR; iR++)
      {
	 recArray[iR] = r1 + iR * dR;
      }
   }
   else
   {
      for (iR = 0; iR < nR; iR++)
      {
	 fscanf(fp, "%f\n", &recArray[iR]);
      }
   }
   fclose(fp);

   /* loop over frequencies */
   do
   {
      masterId = RecvInt(wInfo, 2, -1, FREQUENCY_LIMITS);
      
      if (verbose)
      {
	 fprintf(logInfo.fp_log, 
		 "Frequency limits received: [%d, %d] : [%f, %f]\n", 
		 wInfo[0], wInfo[1], wInfo[0] * dF, wInfo[1] * dF);
	 fflush(logInfo.fp_log);
      }
      
      /* first frequency */
      f1 = wInfo[0] * dF;

      for (indexF = NINT(f1 / dF), f = f1, iF = 0; 
	   iF < (wInfo[1] - wInfo[0] + 1); iF++, f += dF, indexF++)
      {
	 if (verbose)
	 {
	    fprintf(logInfo.fp_log, 
		    "Slave modeling frequency (Hz) : %f: \n", f);
	    fflush(logInfo.fp_log);
	 }
	 
	 /* reseting */
	 for (iR = 0; iR < nR; iR++)
	 {
	    aux11[iR] = zeroC;	         aux12[iR] = zeroC;
	    aux21[iR] = zeroC;	         aux22[iR] = zeroC;
	    aux11Old[iR] = zeroC;	 aux12Old[iR] = zeroC;
	    aux21Old[iR] = zeroC;	 aux22Old[iR] = zeroC;
	    uW[iR][iF] = zeroC;
	 }
	 
	 w = 2 * PI * f;
	 wC.r = w; wC.i = -tau;
	 
	 /* module and phase of complex frequency */
	 wCR = sqrt(wC.r * wC.r + wC.i * wC.i);
	 wCP = atan2(wC.i, wC.r);
	 
	 /* complex slowness step */
	 dUC.r = w * dU / wCR;
	 dUC.i = tau * dU / wCR;
	 
	 /* wCR / wR */
	 wCRwR = wCR / wR;
	 
	 /* auxiliary variable */
	 wCCte.r = wC.r * cte;
	 wCCte.i = wC.i * cte;
	 
	 horSlowness();
	 
	 for (u = u1, iU = 0; iU < nU; iU++, u += dU, 
	      uC.r += dUC.r, uC.i += dUC.i)
	 {
	    uC.r = u;
	    uC.i = u * tau / wRef;
	    
	    uC2.r = 2 * uC.r;
	    uC2.i = 2 * uC.i;
	    
	    aux = uC.r * uC.r - uC.i * uC.i;
	    uuC.i = 2 * uC.r * uC.i;
	    uuC.r = aux;
	    
	    uuC2.r = 2 * uuC.r;
	    uuC2.i = 2 * uuC.i;
	    
	    muC.r = uC.r * -1;
	    muC.i = uC.i * -1;
	    
	    Rm();

	    Rp();
	    
	    /* reseting */
	    As1 = zeroC;      As2 = zeroC;      /* downgoing waves */
	    Cs1 = zeroC;      Cs2 = zeroC;      /* downgoing waves */
	    Bs1 = zeroC;      Bs2 = zeroC;      /* upgoing waves */
	    Ds1 = zeroC;      Ds2 = zeroC;      /* upgoing waves */
	    
	    /* P-wave potential */
	    /* PSlowness^2 - uuC */
	    auxm1 = PSlowness[0].r - uuC.r;
	    auxm2 = PSlowness[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;
	    
	    As1 = uC;
	    if (directWave) Bs1 = muC;
	    
	    /* 1 / am */
	    aux = am.r * am.r + am.i * am.i;
	    amInv.r = am.r / aux;
	    amInv.i = -am.i / aux;
	    
	    /* amInv * uuC */
	    aux2.r = amInv.r * uuC.r - uuC.i * amInv.i;
	    aux2.i = amInv.r * uuC.i + amInv.i * uuC.r;
	    
	    /* aux2 * -I */
	    As2.r = aux2.i;
	    As2.i = -aux2.r;
	    
	    /* notice that Bs2 = As2 */
	    if (directWave) Bs2 = As2;

	    /* S-wave potential */
	    /* SSlowness^2 - uuC */
	    auxm1 = SSlowness[0].r - uuC.r;
	    auxm2 = SSlowness[0].i - uuC.i;
	    
	    /* computing bm */
	    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;
	    
	    /* 1 / bm */
	    aux = bm.r * bm.r + bm.i * bm.i;
	    bmInv.r = bm.r / aux;
	    bmInv.i = -bm.i / aux;
	    
	    /* 1. / bm * uuC */
	    aux1.r = bmInv.r * uuC.r - bmInv.i * uuC.i;
	    aux1.i = bmInv.r * uuC.i + bmInv.i * uuC.r;
	    
	    /* notice that Cs1 = Ds1 */
	    Cs1 = aux1;
	    if (directWave) Ds1 = aux1;
	    
	    Cs2.r = -uC.i;
	    Cs2.i = uC.r;
	    
	    if (directWave)
	    {
	       Ds2.r = -Cs2.r;
	       Ds2.i = -Cs2.i;
	    }
	    
	    /* DD 
	    Bs1 = zeroC;    Bs2 = zeroC;
	    Ds1 = zeroC;    Ds2 = zeroC;*/
	    
	    /* DD 
	    fprintf(logInfo.fp_log, "Ps0.r %f Ps0.i %f u %f f %d \n", 
		    PSlowness[0].r, PSlowness[0].i,
		    uC.r, iF);
	    fprintf(logInfo.fp_log, "As1.r %f As1.i %f\n", As1.r, As1.i);
	    fprintf(logInfo.fp_log, "As2.r %f As2.i %f\n", As2.r, As2.i);
	    fprintf(logInfo.fp_log, "Bs1.r %f Bs1.i %f\n", Bs1.r, Bs1.i);
	    fprintf(logInfo.fp_log, "Bs2.r %f Bs2.i %f\n", Bs2.r, Bs2.i);
	    fprintf(logInfo.fp_log, "Cs1.r %f Cs1.i %f\n", Cs1.r, Cs1.i);
	    fprintf(logInfo.fp_log, "Cs2.r %f Cs2.i %f\n", Cs2.r, Cs2.i);
	    fprintf(logInfo.fp_log, "Ds1.r %f Ds1.i %f\n", Ds1.r, Ds1.i);
	    fprintf(logInfo.fp_log, "Ds2.r %f Ds2.i %f\n", Ds2.r, Ds2.i);*/

	    /* computing compensation for free-surface */
	    /* auxiliar quantities */
	    ambm.r = am.r * bm.r - am.i * bm.i;
	    ambm.i = am.r * bm.i + am.i * bm.r;
	    
	    aux1.r = uC2.r * S2Velocity[0].r - uC2.i * S2Velocity[0].i;
	    aux1.i = uC2.r * S2Velocity[0].i + uC2.i * S2Velocity[0].r;
	    
	    aux2.r = 1. - 
	       (uuC2.r * S2Velocity[0].r - uuC2.i * S2Velocity[0].i);
	    aux2.i = -(uuC2.r * S2Velocity[0].i + uuC2.i * S2Velocity[0].r);
	    
	    /* S2Velocity[0] * S2Velocity[0] */
	    auxm1 = S2Velocity[0].r * S2Velocity[0].r - 
	       S2Velocity[0].i * S2Velocity[0].i;
	    auxm2 = 2 * S2Velocity[0].r * S2Velocity[0].i;
	    
	    /* S2Velocity[0] * S2Velocity[0] * 4 * uuC */
	    aux = 4 * (uuC.r * auxm1 - uuC.i * auxm2);
	    auxm2 = 4 * (uuC.r * auxm2 + uuC.i * auxm1);
	    auxm1 = aux;
	    
	    /* S2Velocity[0] * S2Velocity[0] * 4 * uuC * ambm */
	    aux = ambm.r * auxm1 - ambm.i * auxm2;
	    auxm2 = ambm.r * auxm2 + ambm.i * auxm1;
	    auxm1 = aux;
	    
	    /* aux2 * aux2 */