C++ sf_floatalloc函数代码示例

void fwi(sf_file Fdat, sf_file Finv, sf_file Ferr, sf_file Fgrad, sf_mpi *mpipar, sf_sou soupar, sf_acqui acpar, sf_vec_s array, sf_fwi_s fwipar, sf_optim optpar, bool verb1, int seislet)
/*< fwi >*/
{
	int iter=0, flag;
	float fcost;
	float *x, *direction, *grad;
	sf_gradient gradient;
	FILE *fp;

	/* initialize */
	gradient_init(Fdat, mpipar, soupar, acpar, array, fwipar, verb1);

	/* gradient type */
	gradient=gradient_standard;
	x=array->vv;

	/* calculate first gradient */
	grad=sf_floatalloc(nzx);
	gradient(x, &fcost, grad);

	/* output first gradient */
	if(mpipar->cpuid==0) sf_floatwrite(grad, nzx, Fgrad);

	/* if onlygrad=y, program terminates */
	if(fwipar->onlygrad) return; 

	if(mpipar->cpuid==0) fp=fopen("iterate.txt","a");

	direction=sf_floatalloc(nzx);
	optpar->sk=sf_floatalloc2(nzx, optpar->npair);
	optpar->yk=sf_floatalloc2(nzx, optpar->npair);

	optpar->igrad=1;
	optpar->ipair=0;
	optpar->ils=0;
	optpar->fk=fcost;
	optpar->f0=fcost;
	optpar->alpha=1.;
	/* initialize data error vector */
	for(iter=0; iter<optpar->nerr; iter++){
		optpar->err[iter]=0.;
	}
	optpar->err[0]=optpar->fk;
	if (optpar->err_type==1) optpar->err[optpar->nerr/2]=swap;

	iter=0;
	if(mpipar->cpuid==0){
		l2norm(nzx, grad, &optpar->gk_norm);
		print_iteration(fp, iter, optpar);
	}

	/* optimization loop */
	for(iter=0; iter<optpar->niter; iter++){
		if(mpipar->cpuid==0) sf_warning("-------------------iter=%d----------------------", iter+1);
		
		optpar->ils=0;
		if(iter%optpar->repeat==0) optpar->alpha=1.;

		/* search direction */
		if(iter==0){
			reverse(nzx, grad, direction);
		}else{
			lbfgs_update(nzx, x, grad, optpar->sk, optpar->yk, optpar);
			lbfgs_direction(nzx, grad, direction, optpar->sk, optpar->yk, optpar);
		}

		/* line search */
		lbfgs_save(nzx, x, grad, optpar->sk, optpar->yk, optpar);
		line_search(nzx, x, grad, direction, gradient, optpar, threshold, &flag, mpipar->cpuid, 1);
		optpar->err[iter+1]=optpar->fk;
		if (optpar->err_type==1) optpar->err[optpar->nerr/2+iter+1]=swap;
		
		if(mpipar->cpuid==0){
			l2norm(nzx, grad, &optpar->gk_norm);
			print_iteration(fp, iter+1, optpar);
			fclose(fp); /* get written to disk right away */
			fp=fopen("iterate.txt","a");
		}

		if(mpipar->cpuid==0 && flag==2){
			fprintf(fp, "Line Search Failed\n");
			break;
		}

		if(mpipar->cpuid==0 && optpar->fk/optpar->f0 < optpar->conv_error){
			fprintf(fp, "Convergence Criterion Reached\n");
			break;
		}
	} // end of iter

	if(mpipar->cpuid==0 && iter==optpar->niter){
		fprintf(fp, "Maximum Iteration Number Reached\n");
	}

	/* output vel & misfit */
	if(mpipar->cpuid==0) sf_floatwrite(x, nzx, Finv);
	if(mpipar->cpuid==0) sf_floatwrite(optpar->err, optpar->nerr, Ferr);
	if(mpipar->cpuid==0) fclose(fp);

	return;
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{


  int nt, nx, it, ix, niter, iter, ntfft, nxfft,np, ip, ikt, ikx, iktn, ikxn, ifsnr; /* iktn, ikxn, iNyquist*/
  float dt, dx, pmin, pmax, dp, p, cmax, scalar, sembpmax, num, den;
  float *sembp, *mask, *gy, *fden, *fshift, *SNR;
  float **fdata, **taup, **odata, **tdata, **odatat, **semb; /* tdata is the true data */
  kiss_fft_cpx **cdata, **cdatat;
  char *type;
  sf_file inp, outp, m, spec1, spec2, trued, snr; 


  sf_init(argc,argv);

  inp=sf_input("in");
  m=sf_input("mask");
  outp=sf_output("out");

  if(!sf_histint(inp,"n1",&nt)) sf_warning("No n1 in input");
  if(!sf_histint(inp,"n2",&nx)) sf_warning("No n2 in input");
  if(!sf_histfloat(inp,"d1",&dt)) sf_warning("No n1 in input");
  if(!sf_histfloat(inp,"d2",&dx)) sf_warning("No n2 in input");

  ntfft = 2*kiss_fft_next_fast_size((nt+1)/2);
  nxfft = 2*kiss_fft_next_fast_size((nx+1)/2);
  scalar = 1./(ntfft*nxfft);
  iktn=ntfft/2; ikxn=nxfft/2;
  float dkt = 1.0/(ntfft*dt), fkt = 0.0,kt;
  float dkx = 1.0/(nxfft*dx), fkx = 0.0,kx;

    if (NULL == (type=sf_getstring("type"))) type="amplitude";
    /* [amplitude, semblance] thresholding type, the default is amplitude thresholding  */ 

  	if(!sf_getint("niter",&niter)) niter = 10;
	/* Get the number of iterations */

  	if(!sf_getint("ifsnr",&ifsnr)) ifsnr = 0;
	/* If compute SNR during iteration */


	if(type[0]=='s')
	{

  		if(!sf_getfloat("pmin",&pmin)) pmin=-2;
        /* minimum p */		
  		if(!sf_getfloat("pmax",&pmin)) pmax=2;
        /* maximum p */			
  		if(!sf_getint("np",&np)) np=nx;
        /* number of p */
		dp=(pmax-pmin)/(np-1);		
  		sembp =sf_floatalloc(np);
  		semb  =sf_floatalloc2(nt,np);
 		taup  =sf_floatalloc2(nt,np);	
	}

	/* output files */
  	if (NULL!=sf_getstring("spec2")) 
  	{
		spec2=sf_output("spec2");
		sf_putint(spec2, "n1", ntfft);
		sf_putint(spec2, "n2", nxfft);
	}	       		
  	if (NULL!=sf_getstring("spec1")) 
  	{
		spec1=sf_output("spec1");
		sf_putint(spec1, "n1", ntfft);
		sf_putint(spec1, "n2", nxfft);
	}	  
  	if (ifsnr==1 && (NULL!=sf_getstring("true"))) 
  	{
		snr=sf_output("snr");
		trued=sf_input("true");
		tdata=sf_floatalloc2(nt,nx);
		SNR=sf_floatalloc(niter);
		sf_floatread(tdata[0],nt*nx,trued);

		sf_putint(snr,"n1",niter);
		sf_putint(snr,"d1",1);
		sf_putint(snr,"n2",1);
	}	

  /* Allocate memory */
  cdata =(kiss_fft_cpx**) sf_complexalloc2(ntfft,nxfft);
  cdatat =(kiss_fft_cpx**) sf_complexalloc2(ntfft,nxfft); /* temporary file */
  fshift= sf_floatalloc(ntfft);
  fden 	= sf_floatalloc(ntfft);
  gy 	= sf_floatalloc(nxfft);
  mask  =sf_floatalloc(nx);

  odata =sf_floatalloc2(nt,nx); 
  odatat =sf_floatalloc2(ntfft,nxfft); 
  fdata =sf_floatalloc2(ntfft,nxfft);
  memset(&odata[0][0],0,ntfft*nxfft*sizeof(float)); 

  /* Read data */
  sf_floatread(odata[0],nt*nx,inp);
  sf_floatread(mask,nx,m);

	if(type[0]=='s')
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	sf_oaxa(Fcip,ahx,1);
	sf_oaxa(Fcip,ahy,2);
	sf_oaxa(Fcip,ahz,3);
	sf_oaxa(Fcip,aht,4);
	sf_oaxa(Fcip,ac ,5);
    }

    if (verb){
	sf_raxa(ahx);
	sf_raxa(ahy);
	sf_raxa(ahz);
	sf_raxa(aht);
	sf_raxa(ac);
	sf_raxa(ath);
	sf_raxa(aph);
    }

    if(anis) {
	/* deviation angle */
	if(! sf_getint  ("nps",&nps)) nps=251;
	if(! sf_getfloat("ops",&ops)) ops=-25;
	if(! sf_getfloat("dps",&dps)) dps=0.2;
	aps = sf_maxa(nps,ops,dps);
	sf_setlabel(aps,"ps");
	sf_setunit (aps,"deg");

	if(verb) sf_raxa(aps);
    } else {
	aps = NULL;
    }

    /*------------------------------------------------------------*/
    /* allocate arrays */
    cip = sf_floatalloc4  (sf_n(ahx),sf_n(ahy),sf_n(ahz),sf_n(aht));
    ang = sf_floatalloc2  (sf_n(ath),sf_n(aph));

    /* read velocity */
    vep = sf_floatalloc  (sf_n(ac));    
    sf_floatread(vep,sf_n(ac),Fvel);

    ves = sf_floatalloc  (sf_n(ac));    
    sf_floatread(ves,sf_n(ac),Fvel);
	
    /*------------------------------------------------------------*/
    /* read normals */
    nn  = (vc3d*) sf_alloc(sf_n(ac),sizeof(*nn)); /* normals  */
    vc3dread1(Fnor,nn,sf_n(ac));

    if(anis) {
	/* read anisotropy */
	eps = sf_floatalloc   (sf_n(ac));
	sf_floatread(eps,sf_n(ac),Fani);

	dlt = sf_floatalloc   (sf_n(ac));
	sf_floatread(dlt,sf_n(ac),Fani);

	/* read tilts */	
	tt  = (vc3d*) sf_alloc(sf_n(ac),sizeof(*tt));
	vc3dread1(Ftlt,tt,sf_n(ac));
    }

    /*------------------------------------------------------------*/
    /* in-plane azimuth reference */
    vv.dx=1;
    vv.dy=0;
    vv.dz=0;

int main(int argc, char* argv[])
{
 
    bool hermite_false, hermite_true;
    int n1, n2, npml, pad1, pad2, ns, nw;
    float d1, d2, **v, ds, os, dw, ow;
    sf_complex ****f,  ****obs; 
    sf_file in, out, misfit, source, receiver, record;
    char *order;
    int uts, mts, is, i, j, iw, iter, niter;
    float **recloc;
    float **m_old, **m_new;
    float **d_new, **d_old;
    float **g_old, **g_new; 
    sf_complex ****r_new, ****r_old;
    sf_complex ****Fg; 
    float alpha, beta, gnorm, rnorm; 
    float *datamisfit;

    sf_init(argc, argv);

    in = sf_input("in");
    out = sf_output("out");
    misfit = sf_output("misfit");

    if (!sf_getint("uts",&uts)) uts=0;
//#ifdef _OPENMP
//    mts = omp_get_max_threads();
//#else
    mts = 1;
//#endif
    uts = (uts < 1)? mts: uts;

    hermite_false=false;
    hermite_true=true;
    /* Hermite operator */
    
    if (!sf_getint("npml",&npml)) npml=20;
    /* PML width */
    
    if (!sf_getint("niter",&niter)) niter=0; 
    /* Number of iterations */

    if (NULL == (order = sf_getstring("order"))) order="c";
    /* discretization scheme (default optimal 9-point) */

    fdprep_order(order);

    /* read input dimension */
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input.");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input.");

    if (!sf_histfloat(in,"d1",&d1)) sf_error("No d1= in input.");
    if (!sf_histfloat(in,"d2",&d2)) sf_error("No d2= in input.");

    v = sf_floatalloc2(n1,n2);
    sf_floatread(v[0],n1*n2,in);
    
	/* PML padding */
	pad1 = n1+2*npml;
	pad2 = n2+2*npml;

    /* read receiver */ 
    if (NULL == sf_getstring("receiver")) sf_error("Need receiver="); 
    receiver = sf_input("receiver"); 
    recloc=sf_floatalloc2(n1,n2);
    sf_floatread(recloc[0],n1*n2,receiver);

    /* read source */
    if (NULL == sf_getstring("source")) sf_error("Need source=");
    source = sf_input("source");

    if (!sf_histint(source,"n3",&ns)) sf_error("No ns=.");
    if (!sf_histfloat(source,"d3",&ds)) ds=d2;
    if (!sf_histfloat(source,"o3",&os)) os=0.;

    /* read observed data */
    if (NULL == sf_getstring("record")) sf_error("Need record=");
    record = sf_input("record");

    if (!sf_histint(record,"n4",&nw)) sf_error("No nw=.");
    if (!sf_histfloat(record,"d4",&dw)) sf_error("No dw=.");
    if (!sf_histfloat(record,"o4",&ow)) sf_error("No ow=."); 

    f = sf_complexalloc4(n1,n2,ns,nw);
    obs = sf_complexalloc4(n1,n2,ns,nw);

    sf_complexread(f[0][0][0],n1*n2*ns*nw,source);
    sf_complexread(obs[0][0][0],n1*n2*ns*nw,record);
   
    /* allocate variables */
    m_old = sf_floatalloc2(n1,n2);
    m_new = sf_floatalloc2(n1,n2);
    d_old = sf_floatalloc2(n1,n2);
    d_new = sf_floatalloc2(n1,n2);
    g_old = sf_floatalloc2(n1,n2);
    g_new = sf_floatalloc2(n1,n2);
    
    r_old = sf_complexalloc4(n1,n2,ns,nw);
    r_new = sf_complexalloc4(n1,n2,ns,nw);
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
	int   i;
        float *data;
        char *fn;

        sf_init(argc,argv);

        FILE *Fi;
        sf_file Fo;
        
        cjbsegy *tr;
 
        tr = calloc(sizeof(cjbsegy), 1);

        int nx, ny, nazim, nang, ntau;
        float dazim, dang, dtau, fazim=0.0, fang=0.0, ftau=0;
        int dx, dy, fx=0, fy=0;

        if (!sf_getint("nx",&nx)) nx=101;
        if (!sf_getint("ny",&ny)) ny=101;
        if (!sf_getint("nazim",&nazim)) nazim=8;
        if (!sf_getint("nang",&nang)) nang=21;
        if (!sf_getint("ntau",&ntau)) ntau=101;
        if (!sf_getint("dx",&dx)) dx=1;
        if (!sf_getint("dy",&dy)) dy=1;
        if (!sf_getfloat("dazim",&dazim)) dazim=22.5;
        if (!sf_getfloat("dang",&dang)) dang=2.0;
        if (!sf_getfloat("dtau",&dtau)) dtau=0.002;
        if (!sf_getint("fx",&fx)) fx=1;
        if (!sf_getint("fy",&fy)) fy=1;
        if (!sf_getfloat("ftau",&ftau)) ftau=0;
        if (NULL==(fn=sf_getstring("fn"))) fn="kpstm.ladcig.su.agc";

        /* setup I/O files */
        Fo = sf_output("out");

        if((Fi=fopen(fn,"rb"))==NULL)
        {
           printf("File %s open error!\n",fn);
           exit(0);
        }

        fread(tr,sizeof(cjbsegy),1,Fi);
        int iline0=tr->ep;
        sf_warning("ns=%d dt=%f iLineNo=%d ",tr->ns, tr->dt,iline0);


        if(fseek(Fi, 0L, 2) ==-1)
          printf("input file size unknown; Please specify n2\n");
        int nxy=(int) (ftell(Fi)/((60+ntau)*sizeof(float)));
       
        sf_warning("nxy=%d nx=%d ny=%d ",nxy, nx,ny);
        sf_warning("nazim=%d nang=%d ntau=%d",nazim,nang,ntau);
        sf_warning("dx=%d dy=%d dazim=%f dang=%f dtau=%f",dx,dy,dazim,dang,dtau);
        sf_warning("fx=%d fy=%d fazim=%f fang=%f ftau=%f",fx,fy,fazim,fang,ftau);

        if(nxy!=nx*ny*nazim*nang) {
          sf_warning("nx * ny * nazim * nang != nxy ");
          exit(0);  
         };

        sf_putint(Fo,"n1",ntau);
        sf_putint(Fo,"n2",nang);
        sf_putint(Fo,"n3",nazim);
        sf_putint(Fo,"n4",nx);
        sf_putint(Fo,"n5",ny);
        sf_putfloat(Fo,"d1",dtau);
        sf_putfloat(Fo,"o1",ftau);
        sf_putfloat(Fo,"d2",dang);
        sf_putfloat(Fo,"o2",fang);
        sf_putfloat(Fo,"d3",dazim);
        sf_putfloat(Fo,"o3",fazim);
        sf_putfloat(Fo,"d4",dx);
        sf_putfloat(Fo,"o4",fx);
        sf_putfloat(Fo,"d5",dy);
        sf_putfloat(Fo,"o5",fy);
        sf_putstring(Fo,"label1","z");
        sf_putstring(Fo,"label2","angle");
        sf_putstring(Fo,"label3","azimuth");
        sf_putstring(Fo,"label4","x");
        sf_putstring(Fo,"label5","y");
        sf_putstring(Fo,"unit1","ms");
        sf_putstring(Fo,"unit2","degree");
        sf_putstring(Fo,"unit3","degree");
        sf_putstring(Fo,"unit4","m");
        sf_putstring(Fo,"unit5","m");

        data = sf_floatalloc(ntau);

        rewind(Fi);
        for(i=0;;i++)
        {
          fread(tr,sizeof(cjbsegy),1,Fi);
          if(tr->ep != iline0){
            sf_warning("Read iLineNo=%d finished",iline0);
            iline0=tr->ep;
          }
          fread(data,sizeof(float),ntau,Fi);
          if(feof(Fi))break;
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  int verbose;
  sf_file in=NULL, out=NULL;
  int n1_traces;
  int n1_headers;

  char* header_format=NULL;
  sf_datatype typehead;
  /* kls do I need to add this?  sf_datatype typein; */
  float* fheader=NULL;
  float* intrace=NULL;
  float* fprevheader=NULL;
  int numkeys;
  int ikey;
  char** list_of_keys;
  int *indx_of_keys;
  char* skey;
  int indx_of_skey;
  int skeyvalue;
  bool pkeychanged;
  int itrace=0;

  /*****************************/
  /* initialize verbose switch */
  /*****************************/
  sf_init (argc,argv);

  /* verbose flag controls ammount of print */
  /*( verbose=1 0 terse, 1 informative, 2 chatty, 3 debug ) */
  /* fprintf(stderr,"read verbose switch.  getint reads command line.\n"); */
  if(!sf_getint("verbose",&verbose))verbose=1;
    /* \n
     flag to control amount of print
     0 terse, 1 informative, 2 chatty, 3 debug
  */
  sf_warning("verbose=%d",verbose);
 
  /******************************************/
  /* input and output data are stdin/stdout */
  /******************************************/

  if(verbose>0)fprintf(stderr,"read in file name\n");  
  in = sf_input ("in");

  if(verbose>0)fprintf(stderr,"read out file name\n");
  out = sf_output ("out");

  if (!sf_histint(in,"n1_traces",&n1_traces))
    sf_error("input data not define n1_traces");
  if (!sf_histint(in,"n1_headers",&n1_headers)) 
    sf_error("input data does not define n1_headers");

  header_format=sf_histstring(in,"header_format");
  if(strcmp (header_format,"native_int")==0) typehead=SF_INT;
  else                                       typehead=SF_FLOAT;

  if(verbose>0)fprintf(stderr,"allocate headers.  n1_headers=%d\n",n1_headers);
  fheader = sf_floatalloc(n1_headers);
  fprevheader = sf_floatalloc(n1_headers);
 
  if(verbose>0)fprintf(stderr,"allocate intrace.  n1_traces=%d\n",n1_traces);
  intrace= sf_floatalloc(n1_traces);

  if(verbose>0)fprintf(stderr,"call list of keys\n");
 
  /* this sf_getstring will create parameter descrpiton in the self doc */
  sf_getstring("pkey"); 
  /* \n
     A comma seperated list of primary header keys to monitor to determine 
     gathers.  The trace number in the gather is counted and put in the
     skey header location.
     \n
  */
  list_of_keys=sf_getnstring("pkey",&numkeys);
  /* List of the primary keys monitored to determine gathers. */

  if(list_of_keys==NULL)
    sf_error("The required parameter \"pkey\" was not found.");
  /* I wanted to use sf_getstrings, but it seems to want a colon seperated
     list of keys (eg key=offset:ep:fldr:cdp) and I wanted a comma seperated
     list of keys (eg key=offset:ep:fldr:cdp).
  numkeys=sf_getnumpars("pkey");
  if(numkeys==0)
    sf_error("The required parameter \"pkey\" was not found.");
  fprintf(stderr,"alloc list_of_keys numkeys=%d\n",numkeys);
  list_of_keys=(char**)sf_alloc(numkeys,sizeof(char*)); 
  sf_getstrings("pkey",list_of_keys,numkeys);
  */
  /* print the list of keys */
  if(verbose>1){
    fprintf(stderr,"numkeys=%d\n",numkeys);
    for(ikey=0; ikey<numkeys; ikey++){
      fprintf(stderr,"list_of_keys[%d]=%s\n",ikey,list_of_keys[ikey]);
    }
  }
  if(NULL==(skey=sf_getstring("skey")))
    sf_error("the required parameter \"skey\" was not found");
  /* The name of the secondary key created by the program. */
  
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    bool adj, half, verb, normalize;
    int nt, nx, nh, nh2, ix, ih, iy, i, nn, it, **fold, apt;
    float *trace, **image, **v, rho, **stack, *pp, *off;
    float h, x, t, h0, dh, dx, ti, tx, t0, t1, t2, dt, vi, aal, angle;
    sf_file inp, out, vel, gather, offset;

    sf_init (argc,argv);
    inp = sf_input("in");
    vel = sf_input("vel");  
    out = sf_output("out");
       
    if (!sf_getbool("adj",&adj)) adj=true;
    /* adjoint flag (y for migration, n for modeling) */

    if (!sf_getbool("normalize",&normalize)) normalize=true;
    /* normalize for the fold */
 

    if (!sf_histint(inp,"n1",&nt)) sf_error("No n1=");
    if (!sf_histint(inp,"n2",&nx)) sf_error("No n2=");

    if (!sf_histfloat(inp,"o1",&t0)) sf_error("No o1=");
    if (!sf_histfloat(inp,"d1",&dt)) sf_error("No d1=");
    if (!sf_histfloat(inp,"d2",&dx)) sf_error("No d2=");

    if (adj) {
	if (!sf_histint(inp,"n3",&nh)) sf_error("No n3=");
       
	sf_putint(out,"n3",1);
    } else {
	if (!sf_getint("nh",&nh)) sf_error("Need nh=");
	/* number of offsets (for modeling) */
	
	sf_putint(out,"n3",nh);
    }	

    if (NULL != sf_getstring("gather")) {
	gather = sf_output("gather");
    } else {
	gather = NULL;
    }

    if (!sf_getfloat("antialias",&aal)) aal = 1.0;
    /* antialiasing */

    if (!sf_getint("apt",&apt)) apt = nx;
    /* integral aperture */

    if (!sf_getfloat("angle",&angle)) angle = 90.0;
    /* angle aperture */

    angle = fabsf(tanf(angle*SF_PI/180.0));

    if (!sf_getbool("half",&half)) half = true;
    /* if y, the third axis is half-offset instead of full offset */

    if (!sf_getbool("verb",&verb)) verb = true;
    /* verbosity flag */

    if (!sf_getfloat("rho",&rho)) rho = 1.-1./nt;
    /* Leaky integration constant */

    if (NULL != sf_getstring("offset")) {
	offset = sf_input("offset");
	nh2 = sf_filesize(offset);

	if (nh2 != nh*nx) sf_error("Wrong dimensions in offset");

	off = sf_floatalloc(nh2);	
	sf_floatread (off,nh2,offset);
	sf_fileclose(offset);
    } else {
	if (adj) {
	    if (!sf_histfloat(inp,"o3",&h0)) sf_error("No o3=");
	    if (!sf_histfloat(inp,"d3",&dh)) sf_error("No d3=");
	    sf_putfloat(out,"d3",1.);
	    sf_putfloat(out,"o3",0.);
	} else {
	    if (!sf_getfloat("dh",&dh)) sf_error("Need dh=");
	    /* offset sampling (for modeling) */
	    if (!sf_getfloat("h0",&h0)) sf_error("Need h0=");
	    /* first offset (for modeling) */
	    sf_putfloat(out,"d3",dh);
	    sf_putfloat(out,"o3",h0);
	}
	
	if (!half) dh *= 0.5;

	off = sf_floatalloc(nh*nx);
	for (ix = 0; ix < nx; ix++) {
	    for (ih = 0; ih < nh; ih++) {
		off[ih*nx+ix] = h0 + ih*dh; 
	    }
	}
	offset = NULL;
    }

    v = sf_floatalloc2(nt,nx);
//.........这里部分代码省略.........

int main (int argc, char ** argv)
{
   /* RSF variables */
   bool verb;
   float *dat=NULL, **datf=NULL, **image=NULL, **tr=NULL, **ts=NULL, **trs=NULL;
   float **px=NULL, **pz=NULL;
   sf_file Fin=NULL, Fout=NULL, Ftt=NULL;
   sf_axis a1,a2,a1t,a2t,a3t,ax,az;
   int     nf,ntaper,i1,i2,i2t,itt;
   float   fmax,df,theta,dtheta,tg,tgtap,tmin,xs,xr,eps;
   aalias  aa;
   fslice  tabtt=NULL;
   int     n1,n2,n1t,n2t,n3t,nx,nz;
   float   o1,d1,o2,d2,o1t,d1t,o2t,d2t,o3t,d3t,ox,oz,dx,dz;

   sf_init(argc,argv);

   Fin  = sf_input  ("in" );
   Fout = sf_output ("out" );
   Ftt  = sf_input  ("ttfile" );

   /* axes */
   a1  = sf_iaxa(Fin,1);   /* data */
   a2  = sf_iaxa(Fin,2);   /* data */
   a1t = sf_iaxa(Ftt,1);   /* traveltime */
   a2t = sf_iaxa(Ftt,2);   /* traveltime */
   a3t = sf_iaxa(Ftt,3);   /* traveltime */

   o1  = sf_o(a1);  n1   = sf_n(a1);  d1  = sf_d(a1);
   o2  = sf_o(a2);  n2   = sf_n(a2);  d2  = sf_d(a2);
   o1t = sf_o(a1t); n1t  = sf_n(a1t); d1t = sf_d(a1t);
   o2t = sf_o(a2t); n2t  = sf_n(a2t); d2t = sf_d(a2t);
   o3t = sf_o(a3t); n3t  = sf_n(a3t); d3t = sf_d(a3t);

   /* migration parameters */
   if(! sf_getbool("verb",&verb))     verb = false; /* verbosity flag */
   if(! sf_getfloat("theta",&theta))  theta = 30.;  /* maximum dip */
   if(! sf_getfloat("dtheta",&dtheta))  dtheta = theta/3;  /* taper zone */
   if(dtheta>theta) dtheta=theta;
   if(! sf_getfloat("df",&df))        df = 5.;    /* anti-aliasing sampling */
   if(!sf_getfloat("fmax",&fmax)) fmax=.5/d1;
   if(fmax>(.5/d1)) fmax=.5/d1;
   if (!sf_getint("ntaper",&ntaper)) ntaper=11;
   if(!sf_getfloat("tmin",&tmin)) tmin=3*d1;
   if(!sf_getfloat("xs",&xs)) sf_error("missing xs parameter\n");

   /* image parameters */
   if (!sf_getint("nx",&nx))   nx=n2t;
   if(!sf_getfloat("ox",&ox))  ox=o2t;
   if(!sf_getfloat("dx",&dx))  dx=d2t;
   if (!sf_getint("nz",&nz))   nz=n1t;
   if(!sf_getfloat("oz",&oz))  oz=o1t;
   if(!sf_getfloat("dz",&dz))  dz=d1t;

   /* checking dimensions */
   if((dx!=d2t)||(dz!=d1t)) 
     sf_error("sampling interval have to be the same in:\n"
	      " image and traveltime file\n");
   if(ox<o2t) ox=o2t; 
   if(oz<o1t) oz=o1t;
   if((ox+(nx-1)*dx)>(o2t+(n2t-1)*d2t)) nx=floor(((o2t+(n2t-1)*d2t)-ox)/dx)+1;
   if((oz+(nz-1)*dz)>(o1t+(n1t-1)*d1t)) nz=floor(((o1t+(n1t-1)*d1t)-oz)/dz)+1;

   /* output axis */
   ax = sf_maxa(nx,ox,dx); if(verb) sf_raxa(ax);
   az = sf_maxa(nz,oz,dz); if(verb) sf_raxa(az);
   sf_oaxa(Fout,az,1);
   sf_oaxa(Fout,ax,2);

   /* anti-aliasing */
   /* df = fmax; */
   nf = initAalias(-1,verb,fmax,df,n1,d1,&aa);
   /* fprintf(stderr,"forcing nf=%d df=%f\n",nf,df); */

   /* aperture angle */
   tg    = tan(SF_PI*theta/180);
   tgtap = tan(SF_PI*(theta-dtheta)/180);
   if(verb) sf_warning("tgmax=%f tgtap=%f",tg,tgtap);

   /* allocating */
   dat   = sf_floatalloc(n1);
   image = sf_floatalloc2(nz,nx);
   datf  = sf_floatalloc2 (n1,nf);
   ts    = sf_floatalloc2(n1t,n2t);
   tr    = sf_floatalloc2(n1t,n2t);
   trs   = sf_floatalloc2(n1t,n2t);
   px    = sf_floatalloc2(n1t,n2t);
   pz    = sf_floatalloc2(n1t,n2t);

   if(verb) sf_warning("initializing traveltime loading");
   /* initializing traveltime maps */
   tabtt = fslice_init(n1t*n2t,n3t,sizeof(float));
   fslice_load(Ftt,tabtt,SF_FLOAT);
   if(verb) sf_warning("traveltime loading has finished");

   /* reading the source-slice from traveltime table */
   eps = .01*d2;
   itt = floor((xs+eps-o3t)/d3t);
   fslice_get(tabtt,itt,ts[0]);
   if(verb) sf_warning("traveltime table from source was read");
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    bool sym;
    float o1,d1, o2,d2, tol;
    int nd, n1, n2, n12, niter, rect1, rect2, nw;
    float **xy, *z, *m, *m2, *d;
    sf_file in, out, coord, pattern;

    sf_init(argc,argv);
    in = sf_input("in");
    out = sf_output("out");
    coord = sf_input("coord");

    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    if (!sf_histint(in,"n1",&nd)) sf_error("No n1= in input");

    if (NULL != sf_getstring("pattern")) {
	/* pattern file for output dimensions */
	pattern = sf_input("pattern");
	
	if (!sf_histint(pattern,"n1",&n1)) sf_error("No n1= in pattern");
	if (!sf_histint(pattern,"n2",&n2)) sf_error("No n2= in pattern");
	if (!sf_histfloat(pattern,"d1",&d1)) d1=1.;
	if (!sf_histfloat(pattern,"d2",&d2)) d2=1.;
	if (!sf_histfloat(pattern,"o1",&o1)) o1=0.;
	if (!sf_histfloat(pattern,"o2",&o2)) o2=0.;
	
	sf_fileclose(pattern);
    } else {
	if (!sf_getint("n1",&n1)) sf_error("Need n1=");
	if (!sf_getint("n2",&n2)) sf_error("Need n2=");
	if (!sf_getfloat("d1",&d1)) d1=1.;
	if (!sf_getfloat("d2",&d2)) d2=1.;
	if (!sf_getfloat("o1",&o1)) o1=0.;
	if (!sf_getfloat("o2",&o2)) o2=0.;
    }

    n12 = n1*n2;

    sf_putint(out,"n1",n1);
    sf_putint(out,"n2",n2);
    sf_putfloat(out,"d1",d1);
    sf_putfloat(out,"d2",d2);
    sf_putfloat(out,"o1",o1);
    sf_putfloat(out,"o2",o2);

    if (!sf_getint("niter",&niter)) niter=10;
    /* number of iterations */

    xy = sf_floatalloc2(2,nd);
    sf_floatread(xy[0],nd*2,coord);

    if (!sf_getint("rect1",&rect1)) rect1=1;
    if (!sf_getint("rect2",&rect2)) rect2=1;
    /* smoothing regularization */
    
    if (!sf_getint("nw",&nw)) nw=2;
    /* interpolator size */
    
    if (!sf_getbool("sym",&sym)) sym=false;
    /* if y, use symmetric shaping */
    if (!sf_getfloat("tol",&tol)) tol=1e-3;
    /* tolerance for stopping iteration */
    
    nnshape_init(sym,nd, n1,n2, o1,o2, d1,d2, 
		 rect1,rect2, nw, 2, xy);
    sf_gmres_init(n12,niter); 
 
    z = sf_floatalloc (n12);
    m = sf_floatalloc (n12);
    d = sf_floatalloc(nd);

    sf_floatread(d,nd,in);
 
    /* make right-hand side */
    nnshape_back(d,z);
    nnshape_smooth(z);
    /* invert */
    sf_gmres(z,m,nnshape,NULL,niter,tol,true);
    if (sym) nnshape_smooth(m);
    
    sf_floatwrite (m, n12, out);

    exit(0);
}

int main(int argc, char* argv[])
{
    int i, niter, n1, n2, n12, i3, n3, rect1, rect2, order;
    float *mm, *dd, **pp, lam;
    bool *known;
    sf_file in, out, dip, mask;

    sf_init (argc,argv);
    in = sf_input("in");
    out = sf_output("out");
    dip = sf_input("dip");

    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
    n12 = n1*n2;
    n3 = sf_leftsize(in,2);

    if (!sf_getint("niter",&niter)) niter=100;
    /* number of iterations */

    if (!sf_getint("order",&order)) order=1;
    /* accuracy order */

    pp = sf_floatalloc2(n1,n2);
    mm = sf_floatalloc(n12);
    known = sf_boolalloc(n12);
    
    if (NULL != sf_getstring ("mask")) {
	mask = sf_input("mask");
	dd = sf_floatalloc(n12);
    } else {
	mask = NULL;
	dd = NULL;
    }

    if (!sf_getint("rect1",&rect1)) rect1=3;
    if (!sf_getint("rect2",&rect2)) rect2=3;
    /* smoothing radius */
    
    pwdsl_init(n1,n2,order,rect1,rect2,0.01);
    pwdsl_set(pp);
    sf_mask_init(known);
    
    for (i3=0; i3 < n3; i3++) {
	sf_warning("slice %d of %d",i3+1,n3);

	sf_floatread(mm,n12,in);

	if (NULL != mask) {
	    sf_floatread(dd,n12,mask);
	} else {
	    dd = mm;
	}

	/* figure out scaling and make known data mask */
	lam = 0.;
	for (i=0; i < n12; i++) {
	    if (dd[i] != 0.) {
		known[i] = true;
		lam += 1.;
	    } else {
		known[i] = false;
	    }
	}
	lam = sqrtf(lam/n12);

	/* read dip */
	sf_floatread(pp[0],n12,dip);

	sf_conjgrad_init(n12, n12, n12, n12, lam, 10*FLT_EPSILON, true, true); 
	sf_conjgrad(NULL,sf_mask_lop,pwdsl_lop,dd,mm,mm,niter);
	sf_conjgrad_close();

	sf_floatwrite (mm,n12,out);
    }

    exit(0);
}

int main (int argc, char *argv[])
{
    bool both;
    int ir, nr, n1,n2,n3,n4, m1, m2, m3, n12, n123, nw, nj1, nj2, i3;
    float *u1, *u2, *p;
    sf_file in, out, dip;
    off_t pos=0;
    allpass ap;

    sf_init(argc,argv);
    in = sf_input ("in");
    dip = sf_input ("dip");
    out = sf_output ("out");

    if (SF_FLOAT != sf_gettype(in) ||
            SF_FLOAT != sf_gettype(dip)) sf_error("Need float type");

    if (!sf_histint(in,"n1",&n1)) sf_error("Need n1= in input");
    if (!sf_histint(in,"n2",&n2)) n2=1;
    if (!sf_histint(in,"n3",&n3)) n3=1;
    n12 = n1*n2;
    n123 = n12*n3;
    nr = sf_leftsize(in,3);

    if (!sf_histint(dip,"n1",&m1) || m1 != n1)
        sf_error("Need n1=%d in dip",n1);
    if (1 != n2 && (!sf_histint(dip,"n2",&m2) || m2 != n2))
        sf_error("Need n2=%d in dip",n2);
    if (1 != n3 && (!sf_histint(dip,"n3",&m3) || m3 != n3))
        sf_error("Need n3=%d in dip",n3);

    if (!sf_getbool("both",&both)) both=false;
    /* if y, compute both left and right predictions */

    if (1 == n3) {
        n4=0;
        if (both) sf_putint(out,"n3",2);
    } else {
        if (!sf_getint("n4",&n4)) n4=2;
        /* what to compute in 3-D. 0: in-line, 1: cross-line, 2: both */

        if (n4 > 2) n4=2;
        if (2==n4) sf_putint(out,"n4",both? 4:2);

        if (0 != n4 || both) {
            sf_unpipe(in,(off_t) n123*sizeof(float));
            pos = sf_tell(in);
        }
    }

    if (!sf_getint("order",&nw)) nw=1;
    /* accuracy */
    if (!sf_getint("nj1",&nj1)) nj1=1;
    /* in-line aliasing */
    if (!sf_getint("nj2",&nj2)) nj2=1;
    /* cross-line aliasing */

    for (ir=0; ir < nr; ir++) {

        if (1 != n4) { /* in-line */
            u1 = sf_floatalloc(n12);
            u2 = sf_floatalloc(n12);
            p  = sf_floatalloc(n12);

            for (i3=0; i3 < n3; i3++) {
                /* read data */
                sf_floatread(u1,n12,in);

                /* read t-x dip */
                sf_floatread(p,n12,dip);

                ap = allpass_init (nw,nj1,n1,n2,1,p);

                /* apply */
                allpass1(false, false, ap, u1, u2);

                /* write t-x destruction */
                sf_floatwrite(u2,n12,out);
            }

            free(u1);
            free(u2);
            free(p);

        }

        if (0 != n4) { /* cross-line */
            u1 = sf_floatalloc(n123);
            u2 = sf_floatalloc(n123);
            p  = sf_floatalloc(n123);

            /* read data */
            sf_seek(in,pos,SEEK_SET);
            sf_floatread(u1,n123,in);

            /* read t-y dip */
            sf_floatread(p,n123,dip);

            ap = allpass_init(nw,nj2,n1,n2,n3,p);

//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    int i, iter, niter, p[6][2], status, *mask;
    float *buf, *buf2, *wht;
    double rn, rnp, alpha, beta;
    pid_t pid[6]={1,1,1,1,1,1};
    off_t nm, nd, msiz, dsiz, pos;
    size_t nbuf, mbuf, dbuf;
    FILE *xfile, *Rfile, *gfile, *sfile, *Sfile;
    char *x, *R, *g, *s, *S, *prog;
    sf_file mod, dat, from, mwt, x0, known;  /* input */
    sf_file to, out; /* output */

    extern int fseeko(FILE *stream, off_t offset, int whence);
    extern off_t ftello (FILE *stream);

    sf_init(argc,argv);

    dat = sf_input("in");
    mod = sf_input("mod");

    if (SF_FLOAT != sf_gettype(mod) ||
	SF_FLOAT != sf_gettype(dat)) 
	sf_error("Need float type in mod and dat");

    for (i=0; i < argc-1; i++) {
	argv[i]=argv[i+1];
    }
    for (i=0; i < argc-1; i++) {	
	/* find the program to run */
	if (NULL == strchr(argv[i],'=')) {
	    /* first one without the '=' */
	    prog = argv[0];
	    argv[0] = argv[i];
	    argv[i] = prog;
	    break;
	}
    }

    argv[argc-1] = sf_charalloc(6);
    snprintf(argv[argc-1],6,"adj=X");

    if (!sf_getint("niter",&niter)) niter=1;
    /* number of iterations */

    Rfile = sf_tempfile(&R,"w+b"); 
    xfile = sf_tempfile(&x,"w+b"); 
    gfile = sf_tempfile(&g,"w+b");
    sfile = sf_tempfile(&s,"w+b");
    Sfile = sf_tempfile(&S,"w+b");

    fclose(Rfile);
    fclose(xfile);
    fclose(gfile);
    fclose(sfile);
    fclose(Sfile);

    nm = sf_filesize(mod);
    nd = sf_filesize(dat);

    /* I/O buffers */
    nbuf = BUFSIZ/sizeof(float);
    buf  = sf_floatalloc(nbuf);
    buf2 = sf_floatalloc(nbuf);

    if (NULL != sf_getstring("mwt")) {
	mwt = sf_input("mwt"); /* model weight */
	wht = sf_floatalloc(nbuf);
    } else {
	mwt = NULL;
	wht = NULL;
    }

    if (NULL != sf_getstring("known")) {
	known = sf_input("known"); /* known model mask */
	if (SF_INT != sf_gettype(known)) sf_error("Need int type in known");
	mask = sf_intalloc(nbuf);
    } else {
	known = NULL;
	mask = NULL;
    }

    if (NULL != sf_getstring("x0")) {
	x0 = sf_input("x0"); /* initial model */
    } else {
	x0 = NULL;
    }

    for (i=0; i < 6; i++) { /* make six pipes */
	if (pipe(p[i]) < 0) sf_error("pipe error:");
    }

    for (iter=0; iter < niter; iter++) {
	for (i=0; i < 6; i++) { /* fork six children */
	    if ((pid[i] = fork()) < 0) sf_error("fork error:");
	    if (0 == pid[i]) break;
	}
	
	if (0 == pid[0]) {	
	    /* feeds rr to p[0] */
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
  int nx, ny, nt, npara, nc;
  float x, y, dx, dy, dt, x0, y0, tini, t0;
  float **t0sq,**time, *coeff, ***dtime; 
  float w1,w2,w3,A1,A2,A3,A4,A5,B1,B2,B3,C1,C2,C3,C4,C5,A,B,C;
  int i,j,k,test;
  int count = 0;
  bool verb;

  sf_file inp, out, fit ,inicoef;

  sf_init(argc, argv);
    inp = sf_input("in");
    inicoef = sf_input("coef");
    out = sf_output("out");
    fit = sf_output("fit");

    /* input vector t0^2' */
    if (!sf_histint(inp,"n1",&nx)) sf_error("No n1=");
    if (!sf_histint(inp,"n2",&ny)) sf_error("No n2=");
    if (!sf_histint(inp,"n3",&nt)) sf_error("No n3=");
    if (!sf_histfloat(inp,"d1",&dx)) sf_error("No d1=");
    if (!sf_histfloat(inp,"d2",&dy)) sf_error("No d2=");
    if (!sf_histfloat(inp,"d3",&dt)) sf_error("No d3=");
    if (!sf_histfloat(inp,"o1",&x0)) sf_error("No o1=");
    if (!sf_histfloat(inp,"o2",&y0)) sf_error("No o2=");
    if (!sf_histfloat(inp,"o3",&tini)) sf_error("No o2=");

    /*Number of fitting parameters*/
    npara=16;
    
    /*Verbal flag*/
    if (!sf_getbool("verb",&verb)) verb=false;
    
    /*Memory allocation*/
    coeff = sf_floatalloc(npara);
    t0sq = sf_floatalloc2(nx,ny);
    time = sf_floatalloc2(nx,ny);
    dtime = sf_floatalloc3(nx,ny,npara);
    
    /*Output dimension*/
    if (!sf_histint(inicoef,"n1",&nc) || nc != npara) 
	sf_error("Need n1=%d in inicoef",npara);
    
    /*Shift the third dimension to 4th to insert coefficients*/
     sf_shiftdim(inp, fit, 3);
     sf_putint(fit,"n3",npara);
     sf_putint(fit,"d3",1);
     sf_putint(fit,"o3",0);
    
    /* Loop over time slices */
      for(k=0; k < nt; k++) {
        /*Read intial parameters*/
        sf_floatread(coeff,npara,inicoef);
        w1 = coeff[0];w2 = coeff[1];w3 = coeff[2];
        A1 = coeff[3];A2 = coeff[4];A3 = coeff[5];A4 = coeff[6];A5 = coeff[7];
        B1 = coeff[8];B2 = coeff[9];B3 = coeff[10];
        C1 = coeff[11];C2 = coeff[12];C3 = coeff[13];C4 = coeff[14];C5 = coeff[15];
        sf_floatread(t0sq[0],nx*ny,inp);
        

        /* Loops for each x and y*/
            for(j=0;j<ny;j++){
                y = y0 + j*dy;
                
                for(i=0;i<nx;i++){
                    x = x0 + i*dx;
                    t0 = sqrt(t0sq[j][i]);
                    
                    /* Avoid dividing by zero*/
                    if (x!=0 || y!=0 || t0!=0) {
                    
                    A = A1*pow(x,4) + A2*pow(x,3)*y + A3*pow(x,2)*pow(y,2) + A4*x*pow(y,3) + A5*pow(y,4);
                    B = B1*pow(x,2) + B2*x*y + B3*pow(y,2);
                    C = C1*pow(x,4) + C2*pow(x,3)*y + C3*pow(x,2)*pow(y,2) + C4*x*pow(y,3) + C5*pow(y,4);
                    
                    /* Compute traveltime (t^2)*/
                    time[j][i] = pow(t0,2) + w1*pow(x,2) + w2*x*y + w3*pow(y,2) + A/(pow(t0,2) + B + sqrt(pow(t0,4) + C + 2*pow(t0,2)*B));
                    /* Compute Derivatives */
                    dtime[0][j][i] = pow(x,2); // Wx
                    dtime[1][j][i] = x*y; // Wxy
                    dtime[2][j][i] = pow(y,2); //Wy

                    dtime[3][j][i] = pow(x,4)/(pow(t0,2) + B + sqrt(pow(t0,4) + C + 2*pow(t0,2)*B)); // A1

                    dtime[4][j][i] = (pow(x,3)*y)/(pow(t0,2) + B + sqrt(pow(t0,4) + C + 2*pow(t0,2)*B)); // A2

                    dtime[5][j][i] = (pow(x,2)*pow(y,2))/(pow(t0,2) + B + sqrt(pow(t0,4) + C + 2*pow(t0,2)*B)); // A3

                    dtime[6][j][i] = (x*pow(y,3))/(pow(t0,2) + B + sqrt(pow(t0,4) + C + 2*pow(t0,2)*B)); // A4

                    dtime[7][j][i] = pow(y,4)/(pow(t0,2) + B + sqrt(pow(t0,4) + C + 2*pow(t0,2)*B)); // A5

                    dtime[8][j][i] = (-A*x*x*(1+pow(t0,2)/(sqrt(pow(t0,4) + C + 2*pow(t0,2)*B))))/pow(pow(t0,2) + B + sqrt(pow(t0,4) + C + 2*pow(t0,2)*B),2); // B1

                    dtime[9][j][i] = (-A*x*y*(1+pow(t0,2)/(sqrt(pow(t0,4) + C + 2*pow(t0,2)*B))))/pow(pow(t0,2) + B + sqrt(pow(t0,4) + C + 2*pow(t0,2)*B),2); // B2

                    dtime[10][j][i] = (-A*y*y*(1+pow(t0,2)/(sqrt(pow(t0,4) + C + 2*pow(t0,2)*B))))/pow(pow(t0,2) + B + sqrt(pow(t0,4) + C + 2*pow(t0,2)*B),2); // B3
                    
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    int n1,n2,n3, jump, i1,i2,i3, j;
    float d, *pp=NULL, *zero=NULL, *one=NULL;
    sf_file in=NULL, out=NULL, mask=NULL;

    sf_init(argc,argv);
    in = sf_input("in");
    out = sf_output("out");
    mask = sf_output("mask");

    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
    n3 = sf_leftsize(in,2);

    if (!sf_getint("jump",&jump)) jump=2;
    /* aliasing */

    if (n2 > 1) {
	sf_putint(out,"n2",n2*jump);
	sf_putint(mask,"n2",n2*jump);
	if (sf_histfloat(in,"d2",&d)) {
	    sf_putfloat(out,"d2",d/jump);
	    sf_putfloat(mask,"d2",d/jump);
	}
    }
    if (n3 > 1) {
	sf_putint(out,"n3",n3*jump);
	sf_putint(mask,"n3",n3*jump);
	if (sf_histfloat(in,"d3",&d)) {
	    sf_putfloat(out,"d3",d/jump);
	    sf_putfloat(mask,"d3",d/jump);
	}
    }

    pp = sf_floatalloc(n1);
    zero = sf_floatalloc(n1);
    one = sf_floatalloc(n1);

    for (i1=0; i1 < n1; i1++) {
	zero[i1] = 0.;
	one[i1] = 1.;
    }

    for (i3=0; i3 < n3; i3++) {
	for (i2=0; i2 < n2; i2++) {
	    sf_floatread(pp,n1,in);
	    sf_floatwrite(pp,n1,out);
	    sf_floatwrite(one,n1,mask);
	    if (n2 >1) {
		for (j=0; j < jump-1; j++) {
		    sf_floatwrite(zero,n1,out);
		    sf_floatwrite(zero,n1,mask);
		}
	    }
	}
	if (n3 > 1) {
	    for (i2=0; i2 < n2*jump; i2++) {
		sf_floatwrite(zero,n1,out);
		sf_floatwrite(zero,n1,mask);
	    }
	}
    }

    exit(0);
}

int main (int argc, char* argv[])
{
    int nh, n1,n2, i1,i2, i, n12, niter, dim, n[SF_MAX_DIM], rect[SF_MAX_DIM];
    int shrt, lng;
    float *trace, *hilb, *num, *den, *phase, mean, c;
    char key[6];
    sf_triangle ts, tl;
    sf_file in, out;

    sf_init (argc,argv);
    in = sf_input("in");
    out = sf_output("out");

    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    
    dim = sf_filedims (in,n);
    n1 = n[0];
    n12 = 1;
    for (i=0; i < dim; i++) {
	snprintf(key,6,"rect%d",i+1);
	if (!sf_getint(key,rect+i)) rect[i]=1;
	n12 *= n[i];
    }
    n2 = n12/n1;

    num = sf_floatalloc(n12);
    den = sf_floatalloc(n12);
    phase = sf_floatalloc(n12);

    if (!sf_getint("niter",&niter)) niter=100;
    /* number of iterations */

    if (!sf_getint("order",&nh)) nh=10;
    /* Hilbert transformer order */
    if (!sf_getfloat("ref",&c)) c=1.;
    /* Hilbert transformer reference (0.5 < ref <= 1) */

    sf_hilbert_init(n1, nh, c);

    if (!sf_getint("short",&shrt)) shrt=1;
    /* short smoothing radius */
    if (!sf_getint("long",&lng)) lng=10;
    /* long smoothing radius */

    ts = sf_triangle_init(shrt,n1,false);
    tl = sf_triangle_init(lng,n1,false);

    mean=0.;
    for (i2=0; i2 < n2; i2++) {
	trace = num+n1*i2;
	hilb = den+n1*i2;

	sf_floatread(trace,n1,in);
	sf_hilbert(trace,hilb);

	for (i1=0; i1 < n1; i1++) {
	    trace[i1] = hypotf(trace[i1],hilb[i1]);
	    hilb[i1] = trace[i1];
	}

	sf_smooth2 (ts, 0, 1, false, trace);
	sf_smooth2 (tl, 0, 1, false, hilb);
	
	for (i1=0; i1 < nh; i1++) {
	    trace[i1] = 0.;
	    hilb[i1] = 0.;
	}
	for (i1=nh; i1 < n1-nh; i1++) {
	    mean += hilb[i1]*hilb[i1];
	}
	for (i1=n1-nh; i1 < n1; i1++) {
	    trace[i1] = 0.;
	    hilb[i1] = 0.;
	}
    }
    mean = sqrtf(n12/mean);

    for (i=0; i < n12; i++) {
	num[i] *= mean;
	den[i] *= mean;
    }

    sf_divn_init(dim, n12, n, rect, niter, true);
    sf_divn (num, den, phase);

    sf_floatwrite(phase,n12,out);

    exit(0);
}