Исправить программу для спрямления (гиперболических) кривых

#include <math.h>
 
#include <rsf.h>
 
#include "fint1.h"
 
static float *vel, *par, h;
static bool slow, squared;
 
static float nmo_map(float t, int it) {
    float v;
    
    v = vel[it];
    if (!squared) v *= v;
    v = slow ? h*v : h/v;
    t = t*t + v;
    if (t > 0.) t=sqrtf(t);
    return t;
}
 
static float shifted_nmo_map(float t, int it) {
    float v, q, s;
    
    v = vel[it];
    s = par[it];
 
    if (!squared) v *= v;
    v = slow ? h*v : h/v;
    q = t*t + s*v;
    if (q > 0.) {
    t=t*(1.-1./s)+sqrtf(q)/s;
    return t;
    } else {
    return q;
    }
}
 
static float taner_nmo_map(float t, int it) {
    float v, a;
    
    v = vel[it];
    a = par[it];
 
    if (!squared) v *= v;
    v = slow ? h*v : h/v;
    t = t*t + v/(1+a*v);
    if (t > 0.) t=sqrtf(t);
    return t;
}
 
int main (int argc, char* argv[])
{
    fint1 nmo;
    bool half;
    int ix,ih, nt,nx,nw, nh, noff, nmask, CDPtype, mute, *mask;
    float dt, t0, h0, dh, dy, str;
    float *trace, *off;
    mapfunc nmofunc;
    sf_file cmp, nmod, velocity, offset, msk, het;
 
    sf_init (argc,argv);
    cmp      = sf_input("in");
    velocity = sf_input("velocity");
    nmod     = sf_output("out");
 
    if (NULL != sf_getstring("s")) {
        het = sf_input("s");
    nmofunc = shifted_nmo_map;  
    } else if (NULL != sf_getstring("a")) {
        het = sf_input("a");
        nmofunc = taner_nmo_map;
    } else {
    het = NULL;
    nmofunc =  nmo_map;
    } 
 
    if (SF_FLOAT != sf_gettype(cmp)) sf_error("Need float input");
    if (!sf_histint  (cmp,"n1",&nt)) sf_error("No n1= in input");
    if (!sf_histfloat(cmp,"d1",&dt)) sf_error("No d1= in input");
    if (!sf_histfloat(cmp,"o1",&t0)) sf_error("No o1= in input");
    if (!sf_histint  (cmp,"n2",&nh)) sf_error("No n2= in input");
    nx = sf_leftsize(cmp,2);
 
    if (!sf_getbool("half",&half)) half=true;
    /* if y, the second axis is half-offset instead of full offset */
    if (!sf_getfloat("str",&str)) str=0.5;
    /* maximum stretch allowed */
 
    if (!sf_getint("mute",&mute)) mute=12;
    /* mute zone */
 
    CDPtype=1;
    off = sf_floatalloc(nh);
 
    if (NULL != sf_getstring("offset")) {
    offset = sf_input("offset");
    if (SF_FLOAT != sf_gettype(offset)) sf_error("Need float offset");
    noff = sf_filesize(offset);
    if (noff == nh) {
        sf_floatread (off,nh,offset);
    } else if (noff != nh*nx) {
        sf_error("Wrong dimensions in offset");
    }
    } else {
    if (!sf_histfloat(cmp,"d2",&dh)) sf_error("No d2= in input");
    if (!sf_histfloat(cmp,"o2",&h0)) sf_error("No o2= in input");
 
    if (sf_histfloat(cmp,"d3",&dy) && !sf_getint("CDPtype",&CDPtype)) {
        CDPtype=half? 0.5+dh/dy : 0.5+0.5*dh/dy;
        if (CDPtype < 1) {
        CDPtype=1;
        } else if (1 != CDPtype) {
        sf_histint(cmp,"CDPtype",&CDPtype);
            sf_warning("CDPtype=%d",CDPtype);
        }
    }       
 
    for (ih = 0; ih < nh; ih++) {
        off[ih] = h0 + ih*dh; 
    }
 
    noff = nh;
    offset = NULL;
    }
    
    if (NULL != sf_getstring("mask")) {
    msk = sf_input("mask");
 
    if (SF_INT != sf_gettype(msk)) sf_error("Need integer mask");
    nmask = sf_filesize(msk);
    mask = sf_intalloc(nh);
 
    if (nmask == nh) {
        sf_intread (mask,nh,msk);
    } else if (nmask != nh*nx) {
        sf_error("Wrong dimensions in mask");
    }
    } else {
    nmask = nh;
 
    msk = NULL;
    mask = NULL;
    }
 
    if (!sf_getbool("slowness",&slow)) slow=false;
    /* if y, use slowness instead of velocity */
 
    if (!sf_getbool("squared",&squared)) squared=false;
    /* if y, the slowness or velocity is squared */
 
    if (!sf_getfloat ("h0",&h0)) h0=0.;
    /* reference offset */
    if (half) h0 *= 2.;
    if (!sf_getint("extend",&nw)) nw=4;
    /* trace extension */
 
    trace = sf_floatalloc(nt);
    vel   = sf_floatalloc(nt);
    if (NULL != het) {
    par = sf_floatalloc(nt);
    } else {
    par = NULL;
    }
 
    nmo = fint1_init (nw, nt, mute);
    
    for (ix = 0; ix < nx; ix++) {
    sf_warning("CMP %d of %d;",ix+1,nx);
 
    sf_floatread (vel,nt,velocity); 
    if (NULL != het) sf_floatread(par,nt,het);
    if (NULL != offset && noff != nh) sf_floatread (off,nh,offset);
    if (NULL != msk && nmask != nh) sf_intread (mask,nh,msk);
 
    for (ih = 0; ih < nh; ih++) {
        sf_floatread (trace,nt,cmp);
        
        /* skip dead traces */
        if (NULL != msk && 0==mask[ih]) {
        sf_floatwrite (trace,nt,nmod);
        continue;
        }
        
        fint1_set(nmo,trace);
 
        h = off[ih];
        if (NULL == offset) h += (dh/CDPtype)*(ix%CDPtype); 
        if (half) h *= 2;
        h = h*h - h0*h0;
        
        stretch(nmo,nmofunc,nt,dt,t0,nt,dt,t0,trace,str);
        sf_floatwrite (trace,nt,nmod);
    }
    }
    sf_warning(".");
 
    exit (0);
}
 
/*  $Id$     */