estimates.cc

COOOL：CWP面向对象最优化库（CWP Object Oriented Optimization Library） COOOL是C++类的一个集合

//
// statics.c 
// 
// This code computes the stacking power for a data set given a set of
// surface consistent residual statics parameters
//
//				authors:	Wenceslau Gouveia
//				October 17, 1994
// **** modified by H.L. Deng, 1996
//

#include <su.h>
#include <header.h>
#include <segy.h>
#include <stream.h>
#include "dataset.h"

segy tr;		// SEGY DATA 

extern "C" { int fgettr(FILE *fp, segy *tp); }

void Error(char* s)
{
	fprintf(stderr,"%s\n", s);
	exit(1);
}

main()
{
	char msg[100];			// error message purposes 
	int *SOURCE, *RECEIVER;		// pointers for source and receiver 
					// statics 
	int ***POINTER;			// pointers for data geometry 
	int *FOLDCMP;			// fold of each CMP	
	int ***POINT_CROSS;      	// pointers for cross correlation 
	int NTRACES;			// number of input traces 
	int shift;			// crosscorrelation lag 
	int ipoint, isource, jsource, itrace, icmp;
	int isample, ireceiver, jreceiver, counter;
	int cmp_current, cmp_last;
	int i, j, k;			// general counters  
	float **CROSS;			// crosscorrelation data 
	double stack;			// computed stacking power 
	double static_j, static_i;
	double temp,ds;                 // sampling time interval
	double *vect;			// used to transfer source and 
					// receiver statics to expect 
	double *grad;                   // gradient vector for the output
	FILE *fp;			// input data file 
	FILE *Xfp;			// cross correlation file 

        fp = fopen(DATA_FILE,"r");

        if (fp == NULL)
        {
                sprintf(msg, "Input: can't open %s", DATA_FILE);
                Error(msg);
        }
	else
	{
		fprintf(stderr, "Data file %s read successfully\n", DATA_FILE); 
	}

        Xfp = fopen(CROSS_FILE,"r");

        if (Xfp == NULL)
        {
                sprintf(msg, "Input: can't open %s", CROSS_FILE);
                Error(msg);
        }
	else
        {
                fprintf(stderr, "Cross correlation file %s read successfully\n",
			 CROSS_FILE); 
        }

	// Memory allocation 
	POINTER = alloc3int(2,MAXFOLD,NCMP);
	SOURCE = alloc1int(NSOURCES);
	RECEIVER = alloc1int(NRECEIVERS);
	FOLDCMP = alloc1int(NCMP);
        CROSS = alloc2float(TOTAL_LAG+1, MAXFOLD*(MAXFOLD-1)/2);
        vect = alloc1double(NSOURCES+NRECEIVERS);
	grad = alloc1double(NSOURCES+NRECEIVERS);

	if (POINTER == NULL ||
	    SOURCE == NULL ||
	    RECEIVER == NULL ||
	    FOLDCMP == NULL ||
	    CROSS == NULL ||
	    vect == NULL)
		Error ("Memory allocation problems. Aborting...");

	// Reading data headers to check number of 
	// source and receiver statics 
	NTRACES = 0;
	isource = 0;
	ireceiver = 0;
	icmp = -1;
	cmp_current = 0;
	cmp_last = 0;
	while (fgettr(fp,&tr)) { 	// reading input trace 
		 cmp_current = tr.cdp;
		if (cmp_current != cmp_last)
		{
			if (icmp >= 0) FOLDCMP[icmp] = itrace;
			icmp++;
			itrace = 0;
			cmp_last = cmp_current;
		}	
		 POINTER[icmp][itrace][0] = tr.sx;	
                 POINTER[icmp][itrace][1] = tr.gx;

		 if (icmp == 0 && itrace == 0)
			SOURCE[isource] = POINTER[icmp][itrace][0];

		else if (!(icmp == 0 && itrace == 0) && 
			   SOURCE[isource] != POINTER[icmp][itrace][0])
		{
			i = isource - 1;
			while (i >= 0 && SOURCE[i] != POINTER[icmp][itrace][0])
				i--;
			if (i < 0)
			{
				isource++;
				SOURCE[isource] = POINTER[icmp][itrace][0];
			}
		}
                if (icmp == 0 && itrace == 0)
                        RECEIVER[ireceiver] = POINTER[icmp][itrace][1];

		else if (!(icmp == 0 && itrace == 0) && 
			   RECEIVER[ireceiver] != POINTER[icmp][itrace][1])
		{
			i = ireceiver - 1;
			while (i >= 0 && 
			       RECEIVER[i] != POINTER[icmp][itrace][1])
				i--;
			if (i < 0)
			{
				ireceiver++;
				RECEIVER[ireceiver] = POINTER[icmp][itrace][1];
			}
		}
		itrace++;
		NTRACES++;
	};
	FOLDCMP[icmp] = itrace;	// last CMP 
	ds = 0.001*tr.dt;
	fclose(fp);

	fprintf(stderr, "%d TRACES were read in %d CMPs\n",NTRACES,icmp+1);
	fprintf(stderr, "%d SOURCE STATICS and %d RECEIVER STATICS\n",
		isource+1, ireceiver+1);

	// Pointer manipulation 
	for (icmp = 0; icmp < NCMP; icmp++)
	{
		for (itrace = 0; itrace < FOLDCMP[icmp]; itrace++) 
		{
			isource = 0;
			while (SOURCE[isource] != POINTER[icmp][itrace][0])
                                isource++;

			POINTER[icmp][itrace][0] = isource;

			ireceiver = 0;
			while (RECEIVER[ireceiver] != POINTER[icmp][itrace][1])
				ireceiver++;
				POINTER[icmp][itrace][1] = ireceiver;
		}
	}

	// Freeing some memory 
	free1int(SOURCE);
	free1int(RECEIVER);

	// Defining vector structure for communication with Expect

	vect = alloc1double(NSOURCES+NRECEIVERS);
	int iflag;
	
	rewind(Xfp);
	fread(CROSS[0],sizeof(float),
	      (TOTAL_LAG+1)*(MAXFOLD*(MAXFOLD-1)/2),Xfp);
	fclose(Xfp);

	while (cin >> iflag)
	{
	   for (i=0; i<NSOURCES+NRECEIVERS; i++)
	   {
	      cin >> vect[i];
	      grad[i] = 0.;
	   }

	   for (stack = 0., icmp = 0; icmp < NCMP; icmp++)
	   {
	      for (ipoint = 0, i = 0; 
		   i < FOLDCMP[icmp] - 1; i++)
	      {
		 isource = POINTER[icmp][i][0];
		 ireceiver = POINTER[icmp][i][1];
		 static_i = vect[isource] +  
		    vect[NSOURCES+ireceiver];

		 for (j =i + 1; 
		      j < FOLDCMP[icmp]; j++, ipoint++)
		 {
		    jsource = POINTER[icmp][j][0];
		    jreceiver = POINTER[icmp][j][1];
		    static_j = vect[jsource] + 
		       vect[NSOURCES+jreceiver];

		    shift = TOTAL_LAG/2 + 
		       NINT((static_j - static_i)/ds);

		    if (iflag)
		    {
		       temp = -(double)((CROSS[ipoint][shift+1]
					  - CROSS[ipoint][shift-1])/(2*ds));
		       grad[isource] -= temp;
		       grad[jsource] += temp;
		       grad[NSOURCES+ireceiver] -= temp;
		       grad[NSOURCES+jreceiver] += temp;
		    }
		    else 
		       stack -= (double) CROSS[ipoint][shift];
		     
		    }
		 }
	      }
	      // Exporting the computing stacking power 
	      // to the library via expect 
	      if (iflag)
	      {
		 for (i=0; i<NSOURCES+NRECEIVERS; i++)
		    cout << grad[i] << "  ";
		 cout << endl;
	      }
	      else
	      {
		 cout << stack << endl;
//		 cerr << "\n Stacking Power: " << stack << endl;
	      }
	}
	
     }