view release on metacpan or  search on metacpan

lib/App/SeismicUnixGui/c/synseis/src/synseis.c  view on Meta::CPAN

			source_filename = &argv[i][2];
			resample_source = FALSE;
			if( verbose == TRUE ) {
		          printf("%s\n",source_filename);
			}
		break;

		case 'V': /*-V verbose*/
			  verbose = TRUE;
			  if( verbose == TRUE ) {
		          printf("%s\n",zrhov_filename);
			  }
		break;

		case 'W': /*water depth*/
			sscanf (&argv[i][2],"%f",&water_depth);
			/* not like by linux
			  water_depth=(float)atof(&argv[i][2]);*/
		break;

		case 'X': /* -X0 t0 for resampling source */
			if (argv[i][2]== '0') {
				sscanf (&argv[i][3],"%f",&xstart); /*default=0.0*/
				/* doesn't work in linux*/
				/* 	xstart = (float)atof(&argv[i][3]);defaullt 0 */
			}
		break;

		case 'Z': /*-Zdepth-density-velocityfile name*/
					zrhov_filename = &argv[i][2];
		break;
	   }	/*switch*/
    } /* end '-' */
  } /* end for */

  if( verbose == TRUE ) printf("%s\n",zrhov_filename);
  if( verbose == TRUE ) printf("t_increment= %f\n",t_increment);
  if( verbose == TRUE ) printf("z_increment =%f\n",z_increment);
  if( verbose == TRUE ) printf("water_depth =%f\n",water_depth);
  if( verbose == TRUE ) printf("source_filname =%s\n",source_filename);

/*  if (argc == 1 || error ) { */
/*    fprintf (stderr, " -AFfrequency of Ricker wavelet(Hz)\n"); */
/*    fprintf (stderr, " -AElength of ricker wavelet in secs\n"); */
/*    fprintf (stderr, " -Aoname of output ricker wavelet file if wanted\n");*/
/*    fprintf (stderr, " -CZoutput file with depth, refle. coef. pairs*\n"); */
/*    fprintf (stderr, " -CToutput file with time, reflec. coef. pairs*\n");*/
/*    fprintf (stderr, " -Isampling interval (s)*\n");*/
/*    fprintf (stderr, " -IZsampling interval in depth (m)*\n");*/
/*    fprintf (stderr, " -LDresampled log density filename*\n");*/
/*    fprintf (stderr, " -LVresampled velocity filename*\n");*/
/*    fprintf (stderr, " -Roresampled output source filename*\n");*/
/*    fprintf (stderr, " -Rresample source=TRUE, otherwise don't resample\n");*/
/*    fprintf (stderr, " -Ssource_file_name (resampling=FALSE)\n");*/
/*    fprintf (stderr, " -Vspill all the details of the modelin\n");*/
/*    fprintf (stderr, " -Zname of file with depth(mbsf),density(g/cc),velocity(m/s)*\n");*/
/*    fprintf (stderr, "\n\n\n\n\n * ALWAYS USE\n");*/
/* */
/*  } */

  /*DATA INPUT */

  /*SOURCE IN A FILE */
  /* read source file*/
if(read_source == TRUE) {

	if(  (fpin = fopen(source_filename,"r") ) == NULL) {
		printf("Can't open file1 %s, try again\n", source_filename);
		exit(0);
	}

     /*RESAMPLE SOURCE*/
	if(resample_source == TRUE) {
		for(i=0, num_pts_src=0; (!feof(fpin));  i++, num_pts_src++) {
			fscanf (fpin, "%d %f", &t, &A);
			/*wint[i] = t;*/
			wint[i] = (float)t * t_increment; /*TINT;default 1 ms Sample Interval*/
			winA[i] = A;
		}
		fclose(fpin);
	}

    /* DO NOT RESAMPLE SOURCE */
    if(resample_source == FALSE) {
    	for(i=0, num_pts_src=0; (!feof(fpin));  i++, num_pts_src++) {
    		fscanf (fpin, "%f\n", &A); /* N.B. only amplitudes */
    		winA_reg[i] = A;
    		/*	printf(" amplitude from file= %f\n",winA_reg[i]); */
    	}
    	fclose(fpin);
    	num_pts_reg = num_pts_src;
    	xstart = (float)XSTART; /* XSTART default = 0.0 */
    	/*t_increment = (float)TINT; TINT=0.001 s */

    	if(verbose == TRUE) {
    		printf("For source resampled before i/p \n");
    		printf("t_increment= %f  xstart= %f counter= %f\n",t_increment, xstart,(float)1);

    		for(i=0; i< num_pts_reg; i++) {
    			wint_reg = xstart + (float)i* t_increment;
    			printf("t= %f  A= %f\n",wint_reg, winA_reg[i]);
    		}
    	}
    	/*      t_increment = wint[1] - wint[0]; only amplitudes here
	      if(verbose == TRUE) printf("t_increment=%f\n",t_increment);
	      if(T_INCREMENT != t_increment) printf("t_increment warning\n");*/
	} /* END of DO NOT RESAMPLE SOURCE */

/*resample source */
    if(resample_source == TRUE) {
    	/*regularize source here*/
    	if(verbose == TRUE) printf("    time,      Amplitude\n");
    	subtract = wint[0];
    	for(i=0; i< num_pts_src; i++) {
    		wint[i] = wint[i] - subtract; /* make xstart = 0*/
    		if(verbose == TRUE) printf("%f,      %f\n",wint[i],winA[i]);
    	}
    	regular(xstart, t_increment, num_pts_src,
	      &num_pts_reg, winA, wint, winA_reg); /*xstart from stdin*/

    	if(verbose == TRUE){

lib/App/SeismicUnixGui/c/synseis/src/synseis.c  view on Meta::CPAN

  }

    convolve(Refl_coef_reg_time, num_pts_log_reg, num_pts_reg,
	       t_increment, winA_reg, Convolve_Amplitude);
    num_pts_conv = num_pts_log_reg+num_pts_reg;
    if(verbose == TRUE) {
      for (i=0;  i< num_pts_conv;  i++){
    	  Convolve_time = (t_increment * (float)i  ) + xstart;
    	  printf("%f\t%f\n", Convolve_time, Convolve_Amplitude[i]);
      }
    }
} /* end convolution */
  /*----------------------- my source already resampled ---------------*/

/*read my own source but didn't resample,as it had been already resampled*/
if (read_source == TRUE && resample_source == FALSE){
    /*N.B. num_pts_reg=num_pts_src*/
    xstart = XSTART; /*defaulted to 0.0*/
    /*t_increment = TINT;*/

    if(verbose == TRUE) {
    	printf("own source,didn't resample,already resampled\n");
    	printf("time and Amplitude\n");
    	for(i=0; i < num_pts_src; i++){
    		wint_reg = xstart + t_increment * (float)i;
    		printf("%f      %f\n",wint_reg,winA_reg[i]);
    	}
	 }

    convolve(Refl_coef_reg_time, num_pts_log_reg, num_pts_reg,
	     t_increment, winA_reg, Convolve_Amplitude);
    num_pts_conv = num_pts_log_reg + num_pts_reg;

    if(verbose == TRUE) {
    	printf("time and Amplitude for synthetic seismogram using above source\n");
    	for (i=0;  i< num_pts_conv;  i++){
    		Convolve_time = (t_increment * (float)i  ) + xstart;
    		printf("%f\t%f\n", Convolve_time, Convolve_Amplitude[i]);
      }
    }
  }
  /*---------------  Ricker Source ------------------------------------*/
if (ricker_wavelet == TRUE) { /*calculated ricker source*/
	 t_increment = (float)SI;
	 convolve(Refl_coef_reg_time, num_pts_log_reg, num_samples_Ricker,
	       t_increment, Ricker_Amplitude,Convolve_Amplitude);
      num_pts_conv = num_pts_log_reg + num_samples_Ricker;

/*	if(verbose == TRUE ) {
		printf(" convolving  Ricker Source with reflection coefficient\n");
		for (i=0; i< num_pts_conv;  i++){
			Convolve_time = (t_increment * (float)i  ) + xstart;
			printf("%f\t%f\n", Convolve_time, Convolve_Amplitude[i]);
		}
    } */

 }
  
  /*CONVOLVE SOURCE AND DATA -----------------END-----------------*/
  
  /*OUTPUT CONVOLVED SIGNAL TO STDOUT*/
for (i=0;  i< num_pts_conv;  i++){
	tmin = 2. * water_depth /Vp_water_mps; /*in secs*/
    Convolve_time = (t_increment * (float)i  ) + tmin;
    fprintf(stdout,"%lf\t%lf\n", Convolve_time, Convolve_Amplitude[i] );
 }
	fclose(stdout);

  
} /*end main*/

double ricker(td, freq, rick_long)
  /* subroutine to calculate a ricker wavelet amplitude at a given time value*/
  double 
    td, /*enter in seconds!!!*/
    freq,/* dominant frequency (Hz)*/
    rick_long;   /*length of ricker wavelet*/
  
  {
#define PI        3.1415926
#define PI2       9.8696044
    
    double 
      zero_cross, t2, /* half the signal length*/
      amplitude, /*output value*/
      dominant_freq2; /*dominant frequencey squared, Hz*/
    
    double arg1;
    
    zero_cross = rick_long/2.; /* in seconds*/
    dominant_freq2 = freq * freq; 
    t2 = ( td - zero_cross) *  ( td - zero_cross);
    arg1 = PI2 * dominant_freq2 * t2;
    amplitude = (1. - 2. * arg1) * exp(-arg1);
    return(amplitude);	
  }
  
  
  void regular(xstart,t_increment,num_pts_src,num_pts_reg,winA,wint,winA_reg)
  /*function to regularize data*/

     float xstart, t_increment, winA[], wint[], winA_reg[];
     int num_pts_src,*num_pts_reg;
{	
         /* xstart: x min for resampling
    	t_increment: resampling interval
    	num_pts_src: num_pts in input data set
    	num_pts_reg: output sampled data set
    	winA: input data amplitudes
    	wint: input data times
    	winA_reg: output sampled data amplitudes
    	xint: output interpolation times
    		  first xint has to be >= 0
    		  last xint has to less than the last x value of
    		  the series being resampled*/

  int i; /*local variable*/
  float xint;
  xint =xstart;
/*  printf("made it to subroutine\n"); */
  while(xint < wint[0])  xint = xint +  t_increment;