/***************************************************************/
/*   gsurf.c : generate surface grid from random point data    */
/*              after Smith and Wessel (1990)                  */
/*-------------------------------------------------------------*/
/* Smith, W.H.F. and P. Wessel (1990) Gridding with continuous */
/*  curvature splines in tension. Geophysics, vol.55, 293-305. */
/***************************************************************/
/*  Usage (FORTRAN calling sequence):                          */
/*    call gsurf0(xorg, yorg, smx, smy, mx, my)                */
/*       xorg,yorg : float : coordinates of SW corner of grid  */
/*       smx, smy  : float : mesh size in NS and EW            */
/*       mx, my    : int   : number of grid to N and to E      */
/*           Start and preset fundamental parameters           */
/*    call gsurf1(tens, eps, prelax, maxit)                    */
/*       tens   : float : tension parameter    (default: 0.25) */
/*       eps    : float : convergence limit    (default: 0.)   */
/*       prelax : float : over-relaxation parm (default: 1.4)  */
/*       maxit  : int   : max iteration        (default: 250)  */
/*           Preset optional parameters, if necessary          */
/*    call gsurf2(n, xm, ym, zm)                               */
/*       n      : int         : number of data in the arrays   */
/*       xm, ym : float array : random point coordinates       */
/*       zm     : float array : random point data value        */
/*           Give random data (may be repeated)                */
/*             Calling with n=0 indicates end of random data.  */
/*    call gsurf3(radius, vnul)                                */
/*       radius : float : search radius for trimming           */
/*       vnul   : float : no data value for trimmed-off area   */
/*           Calculate grid value, and trim data               */
/*             If radius <= 0., no trim operation.             */
/*    call gsurf4(fname)                                       */
/*       fname  : char string : filename of grid-data output   */
/*           Output grid data to the file (append)             */
/*             Format is fixed to '((f7.1,9(1x,f7.1)))'        */
/*    call gsurf4m(array, n, m)                                */
/*       array  : float array : area to store grid-data        */
/*       n, m   : int         : adjustable size of array       */
/*           Transfer generated grid data                      */
/*    call gsurf5()                                            */
/*           Close process                                     */
/***************************************************************/
#include    <stdio.h>
#include    <string.h>
#include    <stdlib.h>
#include    <math.h>

static double xorg, yorg, smx, smy, xs, ys;
static double ar, ar2, ar4, tens, tensb, eps, prelax;
static int kseq=0, imx, jmx, imx4, maxit, ntgrd, *kc;
static char *kv; static float *w, *u, *u2, *v, *c;
static double ssr=0., vnul, cts0, cts1, cts2;

int gsurf0(double axo, double ayo, double asmx, double asmy, int mx, int my)
{   int ngrd, n;
    if (kseq != 0)
        { fprintf(stderr, "gsurf0: illegal call sequence\n"); return(-1); }
/* preset fundamental parameters, and prepare work area A (random data) */
    if ((asmx==0.) || (asmy==0.))
        { fprintf(stderr, "gsurf0: mesh size invalid\n"); return(-1); }
    if ((mx<1) || (my<1))
        { fprintf(stderr, "gsurf0: array size rejected\n"); return(-1); }
    xorg = axo; yorg = ayo; imx = mx; jmx = my; imx4 = imx+4;
    smx = asmx; smy = asmy; xs = fabs(smx); ys = fabs(smy);
    ntgrd = imx4*(jmx+4); ngrd = imx*jmx;
    if (ngrd > 10006501)
        { fprintf(stderr, "gsurf0: too huge array\n"); return(-1); }
    if ((w=malloc(ngrd*4*sizeof(float))) == NULL)
        { fprintf(stderr, "gsurf0: memory (w) alloc fail\n"); return(-1); }
    for (n=0; n<(ngrd*4); n+=4) *(w+n) = 1.;
    ar = 1. / (xs/ys); ar2 = ar*ar; ar4 = ar2*ar2;
    tens = 0.25; tensb = tens; eps = 0.; prelax = 1.4; maxit = 250;
    fprintf(stderr, "gsurf0: process start\n");
    kseq = 1; return(0);
}

int gsurf1(double atens, double aeps, double aprelx, int itmax)
{   if (kseq != 1)
        { fprintf(stderr, "gsurf1: illegal call sequence\n"); return(-1); }
/* preset optional parameters */
    if ((atens<0.) || (atens>1.))
        { fprintf(stderr, "gsurf1: tension invalid\n"); return(-1); }
    if (aeps <= 0.)
        { fprintf(stderr, "gsurf1: converge limit invalid\n"); return(-1); }
    if ((aprelx<1.) || (aprelx>=2.))
        { fprintf(stderr, "gsurf1: over-relax parm invalid\n"); return(-1); }
    if (itmax < 1)
        { fprintf(stderr, "gsurf1: max iteration invalid\n"); return(-1); }
    tens = atens; tensb = atens; eps = aeps; prelax = aprelx; maxit = itmax;
    fprintf(stderr, "gsurf1: option paremeters were set\n");
    return(0);
}

int gsurf2(int np, float xp[], float yp[], float zp[])
{   double sx, sy, sxx, sxy, syy, sz, szx, szy, det;
    double fx, fy, zv, dx, dy, rr, xy, xy1, b1, b2, b3, b4, b5;
    float *up, *cp;
    int i, ip, jp, kij, n, ma, mc, m, *kup, *kcp;
    int key;
    if (kseq == 1) kseq = 2;
    if (kseq != 2)
        { fprintf(stderr, "gsurf2: illegal call sequence\n"); return(-1); }
/* read random data, store if usable. (to be repeated) */
    if (np < 0)
        { fprintf(stderr, "gsurf2: data points  n < 0\n"); return(-1); }
    if (np > 0) {
        for (i=0; i<np; i++) {
            fx = (xp[i]-xorg) / smx; fy = (yp[i]-yorg) / smy;
            if ((fx<-0.5) || (fy<-0.5)) continue;
            ip = (int)(fx+0.5); jp = (int)(fy+0.5);
            if ((ip>=imx) || (jp>=jmx)) continue;
            dx = fx - (double)ip; dy = fy - (double)jp; rr = dx*dx + dy*dy;
            kij = (jp*imx + ip) * 4;
            if (rr < *(w+kij)) {
                *(w+kij) = rr; *(w+kij+1) = dx; *(w+kij+2) = dy;
                *(w+kij+3) = zp[i];
            }
        }
        return(0);
    }
    fprintf(stderr, "gsurf2: end of random data, now processing...");
    fflush(stderr);
/* calculate linear trend */
    sx = sy = sxx = sxy = syy = sz = szx = szy = 0.; n = 0;
    for (kij=0,jp=0; jp<jmx; jp++) for (ip=0; ip<imx; ip++,kij+=4) {
        if (*(w+kij) < 1.) {
            fx = *(w+kij+1) + (double)ip; fy = *(w+kij+2) + (double)jp;
            zv = *(w+kij+3); n++; sx += fx; sy += fy;
            sxx += fx*fx; sxy += fx*fy; syy += fy*fy;
            sz += zv; szx += zv*fx; szy += zv*fy;
        }
    }
    det = sxx*syy*n + sx*sy*sxy*2 - sxy*sxy*n - sx*sx*syy - sy*sy*sxx;
    if (det == 0.) { cts0 = cts1 = cts2 = 0.; }
    else {
        cts0 = ( sxx*syy*sz + sx*sxy*szy + sy*sxy*szx
               - sxy*sxy*sz - sx*syy*szx - sy*sxx*szy ) / det;
        cts1 = ( syy*szx*n + sxy*sy*sz + sx*sy*szy
               - sxy*szy*n - sx*syy*sz - sy*sy*szx ) / det;
        cts2 = ( sxx*szy*n + sx*sy*szx + sx*sxy*sz
               - sxy*szx*n - sx*sx*szy - sy*sxx*sz ) / det;
    }
/* remove linear trend */
    mc = 0; ma = 0;
    for (kij=0,jp=0; jp<jmx; jp++) for (ip=0; ip<imx; ip++,kij+=4) {
        if (*(w+kij) < 1.) {
            fx = (double)ip; fy = (double)jp;
            if (*(w+kij) < 0.0025) mc++;
            else { ma++; fx += *(w+kij+1); fy += *(w+kij+2); }
            *(w+kij+3) -= (cts0 + cts1*fx + cts2*fy);
        }
    }
    if ((mc==0) && (ma==0))
        { fprintf(stderr, "\ngsurf2: no data available\n"); return(-1); }
/* rearrange data */
    if ((kc=malloc((mc+1)*sizeof(int))) == NULL)
        { fprintf(stderr, "\ngsurf2: memory (kc) alloc fail\n"); return(-1); }
    if ((c=malloc((mc+1)*sizeof(float))) == NULL)
        { fprintf(stderr, "\ngsurf2: memory (c) alloc fail\n"); return(-1); }
    if ((u=malloc((ma+1)*4*sizeof(float))) == NULL)
        { fprintf(stderr, "\ngsurf2: memory (u) alloc fail\n"); return(-1); }
    kcp = kc; cp = c; up = u;
    for (kij=0,jp=0; jp<jmx; jp++) for (ip=0; ip<imx; ip++,kij+=4) {
        if (*(w+kij) < 1.) {
            if (*(w+kij) < 0.0025) { *kcp++ = kij/4; *cp++ = *(w+kij+3); }
            else {
                kup = (int *)up; *kup = kij/4; up++;
                *up++ = *(w+kij+1); *up++ = *(w+kij+2); *up++ = *(w+kij+3);
            }
        }
    }
/* free work area A, and prepare work area B (grid generation) */
    free(w);
    if ((kv=malloc(ntgrd*sizeof(char))) == NULL)
        { fprintf(stderr, "\ngsurf2: memory (kv) alloc fail\n"); return(-1); }
    for (m=0; m<ntgrd; m++) *(kv+m) = -1;
    if ((v=malloc(ntgrd*sizeof(float))) == NULL)
        { fprintf(stderr, "\ngsurf2: memory (v) alloc fail\n"); return(-1); }
    for (m=0; m<ntgrd; m++) *(v+m) = 0.;
/* prepare for iteration */
    for (m=0; m<mc; m++) {
        kij = *(kc+m); jp = kij/imx; ip = kij%imx;
        kij = (jp+2)*imx4 + (ip+2);
        *(kv+kij) = 5; *(v+kij) = *(c+m);
    }
    free(c); free(kc);
    if ((u2=malloc((ma+1)*2*sizeof(float))) == NULL)
        { fprintf(stderr, "\ngsurf2: memory (u2) alloc fail\n"); return(-1); }
    for (up=u2,m=0; m<ma*4; m+=4) {
        kij = *(int *)(u+m);
        jp = kij/imx; ip = kij%imx; kij = (jp+2)*imx4 + (ip+2);
        dx = *(u+m+1); dy = *(u+m+2); *(v+kij) = *(u+m+3);
        *(kv+kij) = ((dx<0.) ? ((dy<0.)? 3: 4) : ((dy<0.)? 2: 1));
        dx = fabs(dx); dy = fabs(dy); xy = dx + dy; xy1 = xy + 1.;
        b1 = (dy*(xy1+dx) - dx*(1.+dx)*ar2) / (xy*xy1);
        b2 = ((1.+dx)*ar2 - dy) * 2. / xy1;
        b3 = ((1.+dy) - ar2*dx) * 2. / xy1;
        b4 = (dx*(xy1+dy)*ar2 - dy*(1.+dy)) / (xy*xy1);
        b5 = (1.+ar2) * 2. / (xy*xy1);
        *up++ = b1+b2+b3+b4+b5; *up++ = *(u+m+3) * b5;
        *(u+m) = b1; *(u+m+1) = b2; *(u+m+2) = b3; *(u+m+3) = b4;
    }
/* initial grid setting */
    key = 1;
    while (key > 0) {
        key = 0; m = 0;
        for (m=0,jp=0; jp<(jmx+4); jp++) for (ip=0; ip<imx4; ip++,m++) {
            if (*(kv+m) < 0) {
                zv = 0.; rr = 0.;
                if ((jp!=0) && (ip!=0))
                    if (*(kv+m-imx4-1) >= 0)
                        { zv += *(v+m-imx4-1) / 2.; rr += 0.5; }
                if (jp != 0)
                    if (*(kv+m-imx4) >= 0) { zv += *(v+m-imx4); rr += 1.; }
                if ((jp!=0) && (ip!=(imx+3)))
                    if (*(kv+m-imx4+1) >= 0)
                        { zv += *(v+m-imx4+1) / 2.; rr += 0.5; }
                if (ip != 0)
                    if (*(kv+m-1) >= 0) { zv += *(v+m-1); rr += 1.; }
                if (ip != (imx+3))
                    if (*(kv+m+1) >= 0) { zv += *(v+m+1); rr += 1.; }
                if ((jp!=(jmx+3)) && (ip!=0))
                    if (*(kv+m+imx4-1) >= 0)
                        { zv += *(v+m+imx4-1) / 2.; rr += 0.5; }
                if (jp != (jmx+3))
                    if (*(kv+m+imx4) >= 0) { zv += *(v+m+imx4); rr += 1.; }
                if ((jp!=(jmx+3)) && (ip!=(imx+3)))
                    if (*(kv+m+imx4+1) >= 0)
                        { zv += *(v+m+imx4+1) / 2.; rr += 0.5; }
                if (rr == 0.) key++;
                else { *(kv+m) = -2; *(v+m) = zv / rr; }
            }
        }
        for (m=0,jp=0; jp<(jmx+4); jp++) for (ip=0; ip<imx4; ip++,m++)
            { if (*(kv+m) == -2) *(kv+m) = 0; }
    }
    fprintf(stderr, " done\n"); kseq = 3; return(0);
}

int gsurf3(double sr, double avnul)
{   double snet, snetb, denom, cff[12], cfc[12];
    double ccf1, ccf2, c1w0, c1w1, c1s0, c1s1, c2w1, c2w2, c2s1, c2s2;
    double t1, t2, sum, vst, dmax, b1, b2, b3, b4, bsum, bz5;
    float *p, *p2;
    int mx1, mx2, mx3, mx4, mxy, kof[12];
    int i, j, k, ks, n, nit, mit, lit, mi, mj, is, js;
    if (kseq != 3)
        { fprintf(stderr, "gsurf3: illegal call sequence\n"); return(-1); }
/* preset constants for iteration */
    mx1 = imx4; mx2 = mx1+mx1; mx3 = mx2+mx1; mx4 = mx3+mx1;
    mxy = mx1*(jmx+3); snet = 1.-tens; snetb = 1.-tensb;
    denom = (ar4*6.+ar2*8.+6.)*snet + (ar2+1.)*2.*tens;
    cff[4] = cff[7] = -(snet/denom); cff[0] = cff[11] = -(snet/denom)*ar4;
    cff[1] = cff[3] = cff[8] = cff[10] = -(snet/denom)*ar2*2.;
    cff[5] = cff[6] = (((ar2+1.)*snet*4.+tens)/denom);
    cff[2] = cff[9] = (((ar2+1.)*snet*4.+tens)/denom)*ar2;
    cfc[4] = cfc[7] = -snet; cfc[0] = cfc[11] = -snet*ar4;
    cfc[1] = cfc[3] = cfc[8] = cfc[10] = -snet*ar2*2.;
    cfc[5] = cfc[6] = ((ar2+1.)*snet*2.);
    cfc[2] = cfc[9] = ((ar2+1.)*snet*2.)*ar2;
    ccf1 = (1.+ar4)*snet*2.; ccf2 = (1.+ar2)*snet*2. + tens;
    c1w0 = snetb*2. / (snetb+tensb/2.);
    c1w1 = (tensb/2.-snetb) / (snetb+tensb/2.);
    c1s0 = snetb*ar*2. / (snetb*ar + tensb/2.);
    c1s1 = (tensb/2.-snetb*ar) / (snetb*ar + tensb/2.);
    c2w1 = ar2; c2w2 = -(1.+ar2)*2.; c2s1 = 1./ar2; c2s2 = -(1.+1./ar2)*2.;
    kof[0] = +2; kof[1] = +1-mx1; kof[2] = +1; kof[3] = +1+mx1;
    kof[4] = -mx2; kof[5] = -mx1; kof[6] = +mx1; kof[7] = +mx2;
    kof[8] = -1-mx1; kof[9] = -1; kof[10] = -1+mx1; kof[11] = -2;
/* iteration */
    dmax = eps + 1.; nit = 0; lit = 0;
    fprintf(stderr, "gsurf3: iteration - "); fflush(stderr);
    while ((dmax>eps) && (nit<maxit)) {
        p = (v+mx1+2);
        for (i=0; i<imx; i++,p++) { *p = *(p+mx1)*c1w0 + *(p+mx2)*c1w1; }
        p = (v+mxy-mx1+2);
        for (i=0; i<imx; i++,p++) { *p = *(p-mx1)*c1w0 + *(p-mx2)*c1w1; }
        p = (v+mx2+1);
        for (j=0; j<jmx; j++,p+=mx1) { *p = *(p+1)*c1s0 + *(p+2)*c1s1; }
        p = (v+mx3-2);
        for (j=0; j<jmx; j++,p+=mx1) { *p = *(p-1)*c1s0 + *(p-2)*c1s1; }
        p = (v+mx1+1); *p = *(p+2) + *(p+mx2) - *(p+mx2+2);
        p = (v+mx2-2); *p = *(p-2) + *(p+mx2) - *(p+mx2-2);
        p = (v+mxy-mx1+1); *p = *(p+2) + *(p-mx2) - *(p-mx2+2);
        p = (v+mxy-2); *p = *(p-2) + *(p-mx2) - *(p-mx2-2);
        p = (v+2);
        for (i=0; i<imx; i++,p++) { t2 = *(p+mx3) - *(p+mx1);
            t1 = *(p+mx3+1) - *(p+mx1+1) + *(p+mx3-1) - *(p+mx1-1);
            *p = *(p+mx4) + t1*c2w1 + t2*c2w2; }
        p = (v+mxy+2);
        for (i=0; i<imx; i++,p++) { t2 = *(p-mx3) - *(p-mx1);
            t1 = *(p-mx3+1) - *(p-mx1+1) + *(p-mx3-1) - *(p-mx1-1);
            *p = *(p-mx4) + t1*c2w1 + t2*c2w2; }
        p = (v+mx2);
        for (j=0; j<jmx; j++,p+=mx1) { t2 = *(p+3) - *(p+1);
            t1 = *(p+3+mx1) - *(p+1+mx1) + *(p+3-mx1) - *(p+1-mx1);
            *p = *(p+4) + t1*c2s1 + t2*c2s2; }
        p = (v+mx3-1);
        for (j=0; j<jmx; j++,p+=mx1) { t2 = *(p-3) - *(p-1);
            t1 = *(p-3+mx1) - *(p-1+mx1) + *(p-3-mx1) - *(p-1-mx1);
            *p = *(p-4) + t1*c2s1 + t2*c2s2; }
        p=u; p2=u2; dmax = 0.;
        for (k=mx2+2,j=0; j<jmx; j++,k+=4) for (i=0; i<imx; i++,k++) {
            if ((ks=*(kv+k)) == 5) continue;
            if (ks == 0) { sum = 0.;
                for (n=0; n<12; n++) sum += *(v+k+kof[n])*cff[n]; }
            else { t1 = 0.;
                for (n=0; n<12; n++) t1 += *(v+k+kof[n])*cfc[n];
                b1=*p++; b2=*p++; b3=*p++; b4=*p++; bsum=*p2++; bz5=*p2++;
                if (ks == 1) {
                    t2 =  b1*(*(v+k+kof[1])) + b2*(*(v+k+kof[5]))
                        + b3*(*(v+k+kof[9])) + b4*(*(v+k+kof[10])); }
                else if (ks == 2) {
                    t2 =  b1*(*(v+k+kof[3])) + b2*(*(v+k+kof[6]))
                        + b3*(*(v+k+kof[9])) + b4*(*(v+k+kof[8])); }
                else if (ks == 3) {
                    t2 =  b1*(*(v+k+kof[10]))+ b2*(*(v+k+kof[6]))
                        + b3*(*(v+k+kof[2])) + b4*(*(v+k+kof[1])); }
                else if (ks == 4) {
                    t2 =  b1*(*(v+k+kof[8])) + b2*(*(v+k+kof[5]))
                        + b3*(*(v+k+kof[2])) + b4*(*(v+k+kof[3])); }
                sum = (t1 + (t2+bz5)*ccf2) / (ccf1 + ccf2*bsum);
            }
            vst = (sum - *(v+k)) * prelax; *(v+k) += vst;
            if (fabs(vst) > dmax) dmax = fabs(vst);
        }
        nit++; mit = nit*50 / maxit;
        while (lit < mit) { fputc('#', stderr); fflush(stderr); lit++; }
    }
    fprintf(stderr, " done\n"); kseq = 4;
    if (sr <= 0.) return(0);
    vnul = avnul; ssr = sr*sr; mj = (int)(sr/ys);
    fprintf(stderr, "gsurf3: trimming..."); fflush(stderr);
    for (k=mx2+2,j=0; j<jmx; j++,k+=4) for (i=0; i<imx; i++,k++) {
        n = 0;
        for (ks=-mx1,js=-1; (n==0)&&(js<=+mj); js++,ks+=mx1) {
            if (((j+js)<0) || ((j+js)>=jmx)) continue;
            t1 = js*ys; t2 = sqrt(ssr-t1*t1); mi = (int)(t2/xs);
            for (is=-1; is<=+mi; is++) {
                if (((i+is)<0) || ((i+is)>=imx)) continue;
                if (*(kv+k+ks+is) != 0) { n = 1; break; }
            }
        }
        for (ks=+mx1,js=+1; (n==1)&&(js>=-mj); js--,ks-=mx1) {
            if (((j+js)<0) || ((j+js)>=jmx)) continue;
            t1 = js*ys; t2 = sqrt(ssr-t1*t1); mi = (int)(t2/xs);
            for (is=-1; is<=+mi; is++) {
                if (((i+is)<0) || ((i+is)>=imx)) continue;
                if (*(kv+k+ks+is) != 0) { n = 2; break; }
            }
        }
        for (ks=+mx1,js=+1; (n==2)&&(js>=-mj); js--,ks-=mx1) {
            if (((j+js)<0) || ((j+js)>=jmx)) continue;
            t1 = js*ys; t2 = sqrt(ssr-t1*t1); mi = (int)(t2/xs);
            for (is=+1; is>=-mi; is--) {
                if (((i+is)<0) || ((i+is)>=imx)) continue;
                if (*(kv+k+ks+is) != 0) { n = 3; break; }
            }
        }
        for (ks=-mx1,js=-1; (n==3)&&(js<=+mj); js++,ks+=mx1) {
            if (((j+js)<0) || ((j+js)>=jmx)) continue;
            t1 = js*ys; t2 = sqrt(ssr-t1*t1); mi = (int)(t2/xs);
            for (is=+1; is>=-mi; is--) {
                if (((i+is)<0) || ((i+is)>=imx)) continue;
                if (*(kv+k+ks+is) != 0) { n = 4; break; }
            }
        }
        if (n != 4) *(v+k) = vnul;
    }
    fprintf(stderr, " done\n"); return(0);
}

int gsurf4(char *fnam)
{   int i, j, k, n; double t;
    FILE *fp;
    if ((kseq!=4) && (kseq!=5))
        { fprintf(stderr, "gsurf4: illegal call sequence\n"); return(-1); }
/* output generated grid data */
    if (*fnam == '\0') fp = stdout;
    else if ((fp=fopen(fnam,"a")) == NULL) {
        fprintf(stderr, "gsurf4: unable to open file (%s)\n", fnam);
        return (-1);
    }
    k = imx4+imx4+2; n = 0;
    for (j=0; j<jmx; j++,k+=4) {
        for (i=0; i<imx; i++,k++) {
            if ((ssr>0.) && (*(v+k)==vnul)) t = vnul;
            else t = *(v+k) + cts0 + cts1*(double)i + cts2*(double)j;
            if (n != 0) fputc(' ', fp);
            fprintf(fp, "%7.1lf", t);
            if (++n == 10) { fputc('\n',fp); n = 0; }
        }
        if (n > 0) { fputc('\n',fp); n = 0; }
    }
    if (fp != stdout) fclose(fp);
    fprintf(stderr, "gsurf4: grid data were output\n");
    kseq = 5; return(0);
}

int gsurf4m(float *av, int ni, int nj)
{   int i, j, k;
    if ((kseq!=4) && (kseq!=5))
        { fprintf(stderr, "gsurf4m: illegal call sequence\n"); return(-1); }
/* transfer generated grid data */
    k = imx4+imx4+2;
    for (j=0; j<jmx; j++,k+=4) {
        if (j == nj) break;
        for (i=0; i<imx; i++,k++,av++) {
            if (i == ni) break;
            if ((ssr>0.) && (*(v+k)==vnul)) *av = vnul;
            else *av = *(v+k) + cts0 + cts1*(double)i + cts2*(double)j;
        }
        while (i < ni) { av++; i++; }
    }
    fprintf(stderr, "gsurf4m: grid data were tranferred\n");
    kseq = 5; return(0);
}

int gsurf5()
{   if (kseq != 5) {
        fprintf(stderr, "gsurf5: illegal call sequence\n");
        fprintf(stderr, "       forced to set initial state\n");
        kseq = 0; return(1);
    }
    free(u2); free(v); free(kv); free(u);
    fprintf(stderr, "gsurf5: process end\n");
    kseq = 0; return(0);
}

/*******************************************/
/*   FORTRAN Language Interface Routines   */
/*******************************************/
void gsurf0_(float *axo, float *ayo, float *asmx, float *asmy,
             int *mx, int *my)
    { if (gsurf0((double)*axo, (double)*ayo, (double)*asmx, (double)*asmy,
                 *mx, *my) < 0) exit(1); }
void gsurf1_(float *atens, float *aeps, float *aprelx, int *itmax)
    { if (gsurf1((double)*atens, (double)*aeps, (double)*aprelx,
                 *itmax) < 0) exit(1); }
void gsurf2_(int *n, float *xp, float *yp, float *zp)
    { if (gsurf2(*n, xp, yp, zp) < 0) exit(1); }
void gsurf3_(float *sr, float *avnul)
    { if (gsurf3((double)*sr, (double)*avnul) < 0) exit(1); }
void gsurf4_(char *name, size_t lf)
    { char fnam[256]; if (lf > 255) lf=255;
      while ((lf>0) && (*(name+lf-1)==' ')) lf--;
      fnam[0]='\0'; if (lf > 0) strncat(fnam,name,lf);
      if (gsurf4(fnam) < 0) exit(1); }
void gsurf4m_(float *av, int *n, int *m)
    { if (gsurf4m(av, *n, *m) < 0) exit(1); }
void gsurf5_() { gsurf5(); }