view release on metacpan or  search on metacpan

lib/App/SeismicUnixGui/developer/Stripped/model/sufdmod1.su.main.synthetics_waveforms_testpatterns  view on Meta::CPAN

		if (*wfile != '\0' ) warn("Wave file : %s ",wfile);
		else warn("No wave file requested ");
	}

	if ( NINT((float) nz/((float) zd)) + 1  >SU_NFLTS) {
		warn ("Too many depth points : impossible to output wave field. Increase zd ?");
		*wfile='\0';
	}
	

	/* get absorbing conditions
	if (!getparint("abs",abs)) {  abs[0]=0; abs[1]=1;  }
	if (verbose) {
		if (abs[0]==1) warn("absorbing condition on top ");
		if (abs[1]==1) warn("absorbing condition on bottom ");
	}
	/* get decimation coefficients
	if (!getparint("td",&td)) td=1 ;
	if (verbose) warn("time decimation ccoefficent: %d ",td);
	if (!getparint("zd",&zd)) zd=1 ;
	if (verbose) warn("depth decimation ccoefficent: %d ",zd);

	/* choose pressure or particle velocity
	if (!getparint("press", &press)) press=1 ;
	if ((press != 0) && (press != 1)) err ("press must equal 0 or 1");
	if (verbose) {
		if (press==1) warn( "program will output pressure values");
		else if (press==0) warn( "program will output particle velocity values");
	}
		

	/* allocate space
	p=alloc1float(nz);
	v=alloc1float(nz);
	rv=alloc1float(nz);
	rd=alloc1float(nz);
	rd1_5=alloc1float(nz);

	/* read velocity file
	if (*velfile != '\0' ) {
		if ((velocityfp=fopen(velfile,"r"))=='\0') err("cannot open velfile=%s ",velfile);
	}
	if (efread (rv, sizeof(float), nz, velocityfp)!=nz) 
	   err("cannot read %d velocity values ", nz);

	/* read density file  and linearly inderpolate on corrrect location
	if (*dfile != '\0') {
		if ((densityfp=fopen(dfile,"r"))=='\0') err("cannot open dfile=%s ",dfile);
		if (fread(rd,sizeof(float), nz, densityfp)!=nz) err("error reading dfile %s",dfile);
		fclose(densityfp);
	}
	else for (iz=0; iz<nz; iz++) rd[iz]=2500;
	for (iz=0; iz<nz-1; iz++) rd1_5[iz]=(rd[iz]+rd[iz+1])/2;
	rd1_5[nz-1]=rd[nz-1];

	/* time step computation
	vmax=0;
	for (iz=0; iz<nz; iz++) if (rv[iz]>vmax) vmax=rv[iz];if (verbose) warn( "vmax= %f ", vmax);
	dt=dz/1.414/vmax/2; if (verbose) warn( "time step dt= %f ", dt);

	/* maximum number of iterations
	if (nt==0) nt=1+tmax/dt;
	if (verbose) warn( "number of time steps nt= %d ", nt);
	if (NINT( (float) nt/((float)td))+1>SU_NFLTS) err("too many time steps. Increase td ?");

	/* source parameter computation
	   alpha=2*9.8696*freq*freq;

	/* time shift to get a t0 centered source

	if ((styp==0) || (styp == 2)) epst0=fabs(source (0, styp, dt, dz, 0, alpha) / 1e4);
	else if (styp==1) epst0=fabs(source (1/sqrt(2*alpha), styp, dt, dz, 0, alpha)) / 1e4;
	if (verbose) warn( "epst0 = %f ", epst0);

	t=tmax+dt;
	do t=t-dt; while (fabs(source(t, styp, dt, dz, 0, alpha))<epst0);
	t0=t;
	ies=2*t/dt;

	if (verbose) warn("time shift t0 = %f s", t0);

 array initialization
	for (iz=0; iz<nz; iz++) {  v[iz]=0; p[iz]=0;  }

	if (*wfile != '\0') {
		wavefp=fopen (wfile,"w");
		snapsh.d1=dz*zd; snapsh.f1=fz ; snapsh.ns=nz/zd+1; snapsh.d2=dt*td; snapsh.f2=0; 
		/* snapsh.f2=0 is useless since 0 is the "no value" code for SU headers
	}
	/* propagation computation
	itsis=0;
	for (it=0; it<=nt; it++) {
		t=it*dt;
		if (abs[0]==1) p[0]=(p[0]*(1-rv[0]*dt/dz)+2*rd[0]*rv[0]*rv[0]*dt/dz*v[0])/(1+rv[0]*dt/dz);
		else p[0]=0;
		for (iz=1; iz<nz; iz++) p[iz]=p[iz]+rd[iz]*rv[iz]*rv[iz]*dt/dz*(v[iz]-v[iz-1]);
		if (abs[1]!=1) p[nz-1]=0;
		if (it<ies) {
		p[isz]=p[isz]+source(t, styp, dt, dz, t0, alpha);
		}

		for (iz=0; iz<nz-1; iz++) v[iz]=v[iz]+dt/rd1_5[iz]/dz*(p[iz+1]-p[iz]);
		
		if (abs[1] != 1) v[nz-1]=0;
		else
		v[nz-1]=((rd1_5[nz-1]*dz-dt*rd[nz-1]*rv[nz-1])*v[nz-1]-2*dt*p[nz-1])/(rd1_5[nz-1]*dz+dt*rd[nz-1]*rv[nz-1]);

	  if (it % td == 0) {
		   if (press==1) 
			sismo.data[itsis]=p[irz];
		   else
			sismo.data[itsis]=v[irz];
		   itsis++;
		}

		if ((*wfile!='\0') && (it % td == 0)) {
		if (press==1) 
			for (iz=0; iz<nz/zd; ++iz) snapsh.data[iz]=p[iz*zd];
		else
			for (iz=0; iz<nz/zd; ++iz) snapsh.data[iz]=v[iz*zd];