view release on metacpan or  search on metacpan

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

	}
	else {
	  resample_source = TRUE;
/*	  t_increment = (float)atof(&argv[i][2]);*/
	}
	break;
      case 'S': /*-Ssource file name*/
	read_source = TRUE;
	source_filename = &argv[i][2];
	resample_source = FALSE;
	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;
      }
    }
  }
  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 ----------------START-----------------------------*/
//  /* 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);
//    }
//
//    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);
//    }
//
//    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");*/
//    }
//
//    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){
//	printf("  time,     Amplitude,num_pts_reg = %d\n", num_pts_reg);
//      }
//
//      for ( i=0;  i< num_pts_reg;  i++) {
//	wint_reg = xstart + i* t_increment;
//	if(verbose == TRUE) printf("%f,      %f\n",wint_reg, winA_reg[i]);

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

//	printf("%f,%f\n",wint_reg,winA_reg[i]);
//      }
//    }
//*/
//    convolve(Refl_coef_reg, 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]);
//      }
//    }
//  }
//  /*----------------------- 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, 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, 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 /1500.; /*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)
    
     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: in output sampled data set
	winA: input data amplitudes
	wint: input data times
	winA_reg: output sampled data amplitudes
	*/
{	
  /*function to regularize data*/
  int i; /*local variable*/
  float xint;		
  xint =xstart;
  /*printf("made it to subroutine");*/
  while(xint < wint[0])  xint = xint +  t_increment;
  for  (i=0; i < num_pts_src; i++) {
    while(wint[i+1] >= xint) {