posteriori.c

su 的源代码库

   CD = alloc1float(nDM * (nDM + 1) / 2);
   if (PRIOR)
   {
      if(vpFrechet || ipFrechet)
	 CMP = alloc1float(limRange * (limRange + 1) / 2);
      if(vsFrechet || isFrechet)
	 CMS = alloc1float(limRange * (limRange + 1) / 2);
      if(rhoFrechet)
	 CMrho = alloc1float(limRange * (limRange + 1) / 2);
   }
   
   /* FRECHET derivative operator F */
   F = alloc2float(nR * nDM, numberPar * limRange);

   if (IMPEDANCE)
      CmPostInv = 
	 alloc2float(numberParImp * limRange, numberParImp * limRange);
   else
      CmPostInv = alloc2float(numberPar * limRange, numberPar * limRange);

   v1 = alloc2complex(2, numberPar * limRange + 1);
   v2 = alloc2complex(2, numberPar * limRange + 1);
   DmB = alloc3complex(4, numberPar * (limRange + 2), nL);
   derFactor = alloc2complex(2, nL + 1);
   aux11 = alloc2complex(nR, numberPar * limRange);
   aux12 = alloc2complex(nR, numberPar * limRange);
   aux21 = alloc2complex(nR, numberPar * limRange);
   aux22 = alloc2complex(nR, numberPar * limRange);
   aux11Old = alloc2complex(nR, numberPar * limRange);
   aux12Old = alloc2complex(nR, numberPar * limRange);
   aux21Old = alloc2complex(nR, numberPar * limRange);
   aux22Old = alloc2complex(nR, numberPar * limRange);

   /* reading receiver configuration */
   fp = fopen(recFile, "r");
   if (fp == NULL)
   {
      /* standard end-on */
      if (verbose) fprintf(stderr, "No receiver file available\n");
      for (i = 0; i < nR; i++)
      {
         recArray[i] = r1 + i * dR;
      }
   }
   else   
   {
      if (verbose) fprintf(stderr, "Reading receiver file %s\n", recFile);
      for (i = 0; i < nR; i++)
      {
         fscanf(fp, "%f\n", &recArray[i]);
      }
   }
   fclose(fp);
   
   /* reading the model file */
   fp = fopen(modelFile,"r");
   if (verbose)      
     fprintf(stderr,"  Thickness     rho     vP     qP    vS     qS\n");
   for (k = 0; k < nL + 1; k++)
   {
      fscanf(fp, "%f %f %f %f %f %f\n", &thick[k], &rho[k], &alpha[k], 
	     &qP[k], &beta[k], &qS[k]);
      if (verbose)
	fprintf(stderr,"   %7.4f      %4.3f   %3.2f  %5.1f  %3.2f  %5.1f\n",
		thick[k], rho[k], alpha[k], qP[k], beta[k], qS[k]);
   }
   fclose(fp);

   /* setting lim[0] and lim[1] */
   for (depth = thick[0], i = 1; i <= nL; depth += thick[i], i++)
   {
      if (NINT(depth / dZ) <= NINT(limZ[0] / dZ)) lim[0] = i;
      if (NINT(depth / dZ) < NINT(limZ[1] / dZ)) lim[1] = i;
   }
   lim[1]++;

   /* some modeling parameters */
   /* slowness increment */
   dU = (u2 - u1) / (float) nU;

   /* 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 */
   taper = alloc1float(nU);

   /* building window for slowness integration */
   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)));
      }
   }

   /* filtering in frequency domain */
   filter(percW);
   
   /* building frequency filtering */
   /* 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);

   /* correction for the 1st layer */
   thick[0] -= zs;

   /* imaginary part of frequency for damping wrap-around */
   tau = log(tau) / tMod;
   if (tau > TAUMAX)
      tau = TAUMAX;

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

   /* reading data and model covariance matrixes */
   inputCovar(corrDataFile, corrModelFile);
   
   /* starting inverse procedure */
   /* FRECHET derivative matrix  */
      gradient();
   
   if (!noFrechet)
   {
      fp = fopen(frechetFile, "w");
      for (i = 0; i < numberPar * limRange; i++)
      {
	 fwrite(&F[i][0], sizeof(float), nR * nDM, fp);
      }
      fclose(fp);
   }

   /* building a-posteriori model covariance matrix */
   /* prior information is used */
   buffer1 = alloc1float(nDM);
   buffer2 = alloc1float(nDM * nR);

   if (verbose) fprintf(stderr, "Building posteriori covariance...\n");

   for (iParam = 0; iParam < nParam; iParam++)
   {
      for (i = 0; i < nDM; i++)
      {
	 for (offset = i, k = 0; k < nDM; k++)
	 {
	    buffer1[k] = CD[offset];
	    offset += MAX(SGN0(i - k) * (nDM - 1 - k), 1);
	 }

	 /* doing the product CD F */
  	 for (iR = 0; iR < nR; iR++)
	 {
	    buffer2[iR * nDM + i] = 0;
	    for (k = 0; k < nDM; k++)  
	    {
	       buffer2[iR * nDM + i] += buffer1[k] * 
		                        F[iParam][iR * nDM + k];
	    }
 	 }
      }
      
      for (j = 0; j < nParam; j++)
      {
	 CmPostInv[j][iParam] = 0;
	 for (k = 0; k < nDM * nR; k++)
	 {
	    CmPostInv[j][iParam] += buffer2[k] * F[j][k];
	 }
      }
   }

   if (verbose) 
     fprintf(stderr, "Posteriori covariance built. Including prior...\n");

   free1float(buffer1);
   buffer1 = alloc1float(nParam);
   /* including prior covariance matrix */
   if (PRIOR)
   {
       shift = 0;
      if (IMPEDANCE)
      {
	 if (ipFrechet)
	 {
	    for (iParam = 0; iParam < limRange; iParam++)
	    {
	       for (offset = iParam, k = 0; k < limRange; k++)
	       {
		  buffer1[k] = CMP[offset];
		  offset += MAX(SGN0(iParam - k) * (limRange - 1 - k), 1);
	       }
	       
	       for (k = 0; k < limRange; k++) 
	       {
		  CmPostInv[iParam][k] += buffer1[k];
	       }
	    }
            shift += limRange;
	 }
      }
      else
      {
	 if (vpFrechet)
	 {
	    for (iParam = 0; iParam < limRange; iParam++)
	    {
	       for (offset = iParam, k = 0; k < limRange; k++)
	       {
		  buffer1[k] = CMP[offset];
		  offset += MAX(SGN0(iParam - k) * (limRange - 1 - k), 1);
	       }
	       
	       for (k = 0; k < limRange; k++) 
	       {
		  CmPostInv[iParam][k] += buffer1[k];
	       }
	    }
            shift += limRange;
	 }
      }

      if (IMPEDANCE)
      {
	 if (isFrechet)
	 {
	    for (iParam = 0; iParam < limRange; iParam++)
	    {
	       for (offset = iParam, k = 0; k < limRange; k++)
	       {
		  buffer1[k] = CMS[offset];
		  offset += MAX(SGN0(iParam - k) * (limRange - 1 - k), 1);
	       }
	       
	       for (k = 0; k < limRange; k++)
	       {
		  CmPostInv[iParam + shift][k + shift] += buffer1[k];
	       }
	    }
            shift += limRange;
	 }
      }
      else
      {
	 if (vsFrechet)
	 {
	    for (iParam = 0; iParam < limRange; iParam++)
	    {
	       for (offset = iParam, k = 0; k < limRange; k++)
	       {
		  buffer1[k] = CMS[offset];
		  offset += MAX(SGN0(iParam - k) * (limRange - 1 - k), 1);
	       }
	       
	       for (k = 0; k < limRange; k++)
	       {
		  CmPostInv[iParam + shift][k + shift] += buffer1[k];
	       }
	    }
            shift += limRange;
	 }
      }

      if (rhoFrechet)
      {
	 for (iParam = 0; iParam < limRange; iParam++)
	 {
	    for (offset = iParam, k = 0; k < limRange; k++)
	    {
	       buffer1[k] = CMrho[offset];
	       offset += MAX(SGN0(iParam - k) * (limRange - 1 - k), 1);
	    }
	    
	    for (k = 0; k < limRange; k++) 
	    {
	       CmPostInv[iParam + shift][k + shift] += buffer1[k];
	    }
	 }
      }
   }

   if (verbose) fprintf(stderr, "Prior included. Inverting matrix...\n");

   /* freeing memory */
   free1float(buffer1);
   free1float(buffer2);
   free1float(alpha);
   free1float(beta);
   free1float(rho);
   free1float(qP);
   free1float(qS);
   free1float(thick);
   free2complex(PSlowness);
   free2complex(SSlowness);
   free2complex(S2Velocity);
   free1float(CD);
   free1float(CMP);
   free1float(CMS);
   free1float(CMrho);
   free2float(F);
   free2complex(v1);
   free2complex(v2);
   free3complex(DmB);
   free2complex(derFactor); 
   free2complex(aux11);
   free2complex(aux12);
   free2complex(aux21);
   free2complex(aux22); 
   free2complex(aux11Old);
   free2complex(aux12Old);
   free2complex(aux21Old);
   free2complex(aux22Old); 

   /* inverting the matrix */
   CmPost = alloc2float(nParam, nParam);
   for (i = 0; i < nParam; i++) for (j = 0; j < nParam; j++)
      CmPostInv[i][j] = CmPost[i][j];
   inverse_matrix(nParam, CmPostInv);

   if (verbose) fprintf(stderr, "Done with inverse matrix routine.\n");
   
   buffer1 = alloc1float(nParam);
   gp = fopen(postFile, "w");
   for (i = 0; i < nParam; i++)
   {
      fwrite(CmPostInv[i], sizeof(float), nParam, gp);
   }
   fclose(fp);
}