supressure.c

seismic software,very useful

      }

      if( obfp ){
         memcpy( &ob  ,&vi ,240 );
      }

      if( ppfp ){
         memcpy( &pp  ,&vi ,240 );
      }

      if( nhfp ){
         memcpy( &nh  ,&vi ,240 );
      }

      if( chfp ){
         memcpy( &ch  ,&vi ,240 );
      }

      if( dvfp ){
         memcpy( &dv  ,&vi ,240 );
      }

      if( pepfp ){
         memcpy( &pep ,&vi ,240 );
      }

      if( popfp ){
         memcpy( &pop ,&vi ,240 );
      }

      if( mwwfp ){
         memcpy( &mww ,&vi ,240 );
      }

      WD = vi.swdep;

      Z = 0.0;

      j++;

      /*---------------*/
      /* process input */
      /*---------------*/

      for( i=0; i<vi.ns; i++ ){

         if( i == dbg ){
            printf( "break\n" );
         }

         Vobs = vi.data[i];

         if( time ){
            dZout = 0.5 * Vobs * dTout;
         }

         if( Z < Zmin || Z < WD ){
            Vnct = Vobs;
            
         }else{
#if 0
            Vnct = (Vmin*Vmax)/((Vmax-Vmin)*exp(-Kvp*(Z-WD)) + Vmin);
#else
            Vnct = Vmax - (Vmax - Vmin)*exp(-Kvp*(Z-WD));
#endif

         }

         Vp = Vobs *0.3048e-3;
         Vs = A + B *Vp;

         sigma = (Vp*Vp - 2*Vs*Vs) / (2*Vp*Vp - 2*Vs*Vs);

         /*---------------------------------------*/
         /* calculate normal hydrostatic gradient */
         /*---------------------------------------*/

         if( HG ){

            nh.data[i] = HG*Z;
  
         }else{

            nh.data[i] = (P_a0+P_b0*T+P_c0*S)*Z*Z*.3048*.3048 
                + (P_a1+P_b1*T+P_c1*S)*Z*.3048 
                + P_a2;

            /*--------------------------*/
            /* convert from MPa to psia */
            /*--------------------------*/

            nh.data[i] *= 145.0378;

         }

         /*-----------------------------------------*/
         /* compute vertical stress from overburden */
         /*-----------------------------------------*/

         if( Z <= WD ){

            ob.data[i] = nh.data[i];

         }else{

            if( LG ){

               /*---------------------------------*/
               /* compute using constant gradient */
               /*---------------------------------*/

               ob.data[i] = WD*HG + LG*(Z-WD);

            }else{

               /*----------------------------------------*/
               /* compute using integrated density trend */
               /*----------------------------------------*/

               if( Vobs >= Vsalt && Vobs > Vnct ){
                  Rho = 2.12;

               }else{
#if 0
                  Rho = Rho_max - (Rho_max-Rho_min)*exp(-Krho*(Z-WD));
#else
                  Rho = (Rho_max*Rho_min)
                      / ( (Rho_max-Rho_min)*exp(-Krho*(Z-WD))+Rho_min );
#endif

               }

               ob.data[i] = ob.data[i-1] + 0.43353 * Rho * dZout;

            }

         }

         /*---------------------------------*/
         /* account for water depth effects */
         /*---------------------------------*/

         if( Z < WD ){
            c = 0.0;
            Vnct = Vmin;

         }else{
            c = (ob.data[i] - nh.data[i]);

         }

         dv.data[i] = Vnct - Vobs;

         /*------------------------------*/
         /* calculate effective pressure */
         /*------------------------------*/

         if( method == 1 ){

            /*---------------------------*/
            /* normalized Eaton's method */
            /*---------------------------*/

            a = Vobs - Vmin;
            b = Vnct - Vmin;

            if( Vnct > Vmin && Vobs > Vmin ){
               ep.data[i] = c*pow( a / b ,I );

            }else{
               ep.data[i] = 0.0;

            }

         }else{

            /*----------------*/
            /* Eaton's method */
            /*----------------*/

            ep.data[i] = c*pow( Vobs / Vnct ,I );

         }

         /*--------------------------*/
         /* bound effective pressure */
         /*--------------------------*/

         ep.data[i] = ep.data[i] > c ? c : ep.data[i];

         /*-------------------------*/
         /* calculate over pressure */
         /*-------------------------*/

         if( Vobs > Vmin ){
            op.data[i] = c - ep.data[i];

         }else{
            op.data[i] = 0.0;

         }

         /*-------------------------*/
         /* calculate pore pressure */
         /*-------------------------*/

         if( Vobs > Vmin ){


            pp.data[i] = nh.data[i] + op.data[i];
            

         }else{
            pp.data[i] = 0.0;

         }

         if( pp.data[i] > ob.data[i] ){
            pp.data[i] = ob.data[i];
         }

         /*--------------------------------*/
         /* calculate pressure percentages */
         /*--------------------------------*/

         if( (ep.data[i] + op.data[i] ) != 0.0 ){

            pep.data[i] = (100.0 * ep.data[i]) 
                       / (ep.data[i] + op.data[i] );

            pop.data[i] = (100.0 * op.data[i]) 
                       / (ep.data[i] + op.data[i] );

         }else{
            pep.data[i] = 100.0;
            pop.data[i] =   0.0;

         }

/*--------------------------------------------------------------------*\
   calculate fracture pressure based on various models.
\*--------------------------------------------------------------------*/

         if( Vobs > Vmin && Z > WD && Vobs < Vsalt ){

            if( frac == 0 ){

               /* Hubbert-Willis */

               fp.data[i] = ep.data[i]*sigma/(1-sigma) + pp.data[i];

            }else if( frac == 1){

               /* Terzaghi */

               fp.data[i] = 2*ep.data[i]*sigma/(1-sigma) + pp.data[i];

            }else if( frac == 2){

               /* Matthews-Kelly-Constant-Bourgoyne */

               Rms = 1.0 - MKCB_a*exp(-MKCB_b*(Z-WD));
               fp.data[i] = Rms*ep.data[i] + pp.data[i];

            }else if( frac == 3){
#if 0
               /* Fore-Beardsley */
               fp.data[i] = ob.data[i]*(FBmin*FBmax) 
                           / ((FBmax-FBmin)*exp(-FBk*(Z-WD))+FBmin);
#else
               fp.data[i] = nh.data[i] 
                          + c*(FBmax-(FBmax-FBmin)*exp(-FBk*(Z-WD)));
#endif
            }else{

               /* default to no answer*/
               fp.data[i] = ob.data[i];

            }


         }

         if( Vobs >= Vsalt ){

            fp.data[i]  = 0.0;
            pp.data[i]  = 0.0;
            op.data[i]  = 0.0;
            ch.data[i]  = 0.0;

            pop.data[i] = 0.0;
            pep.data[i] = 100.0;

            ep.data[i]  = ob.data[i] - nh.data[i];

         }else if( Vobs < Vsalt && fp.data[i] > 0.0 ){
            ch.data[i] = (fp.data[i] - pp.data[i]) 
                        / (0.43353 * (Rho - Rho_gas ));

         }else{
            ch.data[i] = 0.0;

         }

         if( ch.data[i] > 8000.0 ){
            ch.data[i] = 8000.0;
         }

         mww.data[i] = fp.data[i] - pp.data[i];

         if( ep.data[i]  != ep.data[i]
          || fp.data[i]  != fp.data[i]
          || op.data[i]  != op.data[i]
          || ob.data[i]  != ob.data[i]
          || pp.data[i]  != pp.data[i]
          || nh.data[i]  != nh.data[i]
          || ch.data[i]  != ch.data[i]
          || dv.data[i]  != dv.data[i] 
          || mww.data[i] != mww.data[i]
          || pep.data[i] != pep.data[i]
          || pop.data[i] != pop.data[i] ){
            fprintf( stderr ,"NaN trace=%d sample=%d\n" ,j ,i );
         }

         Z += dZout;
      }

/*--------------------------------------------------------------------*\
            Convert from PSI to PPG if requested
\*--------------------------------------------------------------------*/

      if( ppg ){

         Z = 0;

         for( i=1; i<vi.ns; i++ ){

            if( time ){
               Z += vi.data[i] * 0.5 * dTout;

            }else{
               Z += dZout;
            }

            nh.data[i]  /= Z*0.052;
            ep.data[i]  /= Z*0.052;
            op.data[i]  /= Z*0.052;
            pp.data[i]  /= Z*0.052;
            ob.data[i]  /= Z*0.052;
            fp.data[i]  /= Z*0.052;
            mww.data[i] /= Z*0.052;

         }

      }

      /*----------------*/
      /* output results */
      /*----------------*/

      if( epfp ){
         fputtr(epfp, &ep);
      }

      if( fpfp ){
         fputtr(fpfp, &fp);
      }

      if( pepfp ){
         fputtr(pepfp, &pep);
      }

      if( popfp ){
         fputtr(popfp, &pop);
      }

      if( ppfp ){
         fputtr(ppfp, &pp);
      }

      if( opfp ){
         fputtr(opfp, &op);
      }

      if( obfp ){
         fputtr(obfp, &ob);
      }

      if( dvfp ){
         fputtr(dvfp, &dv);
      }

      if( nhfp ){
         fputtr(nhfp, &nh);
      }

      if( chfp ){
         fputtr(chfp, &ch);
      }

      if( mwwfp ){
         fputtr(mwwfp, &mww);
      }

   } while (fgettr(vifp, &vi));

   return 0;

}

void interpolate( char   method 
                 ,float* xin 
                 ,float* yin 
                 ,float* xout
                 ,float* yout 
                 ,int    nin
                 ,int    nout
                ){

   float yind[4096][4];

   /*----------------------*/
   /* linear interpolation */
   /*----------------------*/

   if( method == 'l' ){
      intlin(nin,xin,yin,yin[0],yin[nin-1],nout,xout,yout);

   /*-------------------------*/
   /* monotonic interpolation */
   /*-------------------------*/

   }else if( method == 'm' ){
      cmonot(nin,xin,yin,yind);
      intcub(0,nin,xin,yind,nout,xout,yout);

   /*---------------------*/
   /* Akima interpolation */
   /*---------------------*/

   }else if( method == 'a' ){
      cakima(nin,xin,yin,yind);
      intcub(0,nin,xin,yind,nout,xout,yout);

   /*----------------------------*/
   /* cubic spline interpolation */
   /*----------------------------*/

   }else if( method == 's' ){
      csplin(nin,xin,yin,yind);
      intcub(0,nin,xin,yind,nout,xout,yout);

   /*--------------------------*/
   /* unknown method specified */
   /*--------------------------*/

   }else{
      err("%s is an unknown interpolation method!\n",method);
   }

}