1 #include <omp.h>
    2 #include <stdio.h>
    3 #include <stdlib.h>
    4 #include <f2c.h>
    5 #include <string.h>
    6 
    7 runloop(p, ij, ip, bnd2, itom, nxmax, nzmax, maxs, maxr, direct,
    8   sh_s0, sh_s0b, sh_s1, sh_rd_glb, sh_array_glb)
    9   integer p;
   10   integer ij, ip, itom, nxmax, nzmax, maxs, maxr;
   11   integer direct;
   12   double bnd2;
   13   sh_ptr_real sh_s0, sh_s0b, sh_s1;
   14   sh_ptr_rd_glb sh_rd_glb;
   15   sh_ptr_array_glb sh_array_glb;
   16 {
   17   integer ig, ii, ix, iz, i_1, i_2, i_3, i_4, i_5;
   18   integer t0m_dim1, xsm_dim1, zsm_dim1, xrm_dim1, zrm_dim1;
   19   integer count_dim1, image_dim1;
   20   integer ir, ns, is, istat, i, rp;
   21   integer nray, icol, kk, k, ncmax, nsamp, nstep;
   22   real xrt, zrt, xst, zst, trad, xsi, zsi, rli, dl1, tci;
   23   real xwavg, xwmin, xwmax, zwmin, zwmax;
   24   real xg, zg, resid;
   25   real thmax, thmin, pi180, pi;
   26   real zsimax, zsimin, tmi, dta, accel;
   27   real error = 0.0;
   28   real r_1, r_2;
   29   ptr_real xray1, zray1, ray1, time, rayl, xray, zray;
   30   real *xw, *zw, *tm;
   31   rd_ptr_real t0m, xsm, zsm, xrm, zrm;
   32   rd_ptr_real s0, s0b, s1;
   33   integer ntot;
   34   rd_ptr_rd_glb rd_glb = sh_rd_glb;
   35   rd_ptr_array_glb array_glb = sh_array_glb;
   36   rd_wr_ptr_real count, image;
   37   rd_ptr_real rc, ri;
   38   integer cur, cur_start, cur_stop;
   39 
   40   extern doublereal xnterp_();
   41 
   42   double atan();
   43   double cos();
   44 
   45   xray1 = (ptr_real) malloc(NMAX * 7 * sizeof(real));
   46   zray1 = &xray1[NMAX];
   47   ray1 = &xray1[2 * NMAX];
   48   xray = &xray1[3 * NMAX];
   49   zray = &xray1[4 * NMAX];
   50   rayl = &xray1[5 * NMAX];
   51   time = &xray1[6 * NMAX];
   52 
   53   s0b = sh_s0b;
   54   if (ip == 1) {
   55     s0 = sh_s0;
   56   } else {
   57     s1 = sh_s1;
   58   }
   59   xw = (real *) malloc(maxs * sizeof(real));
   60   zw = (real *) malloc(maxs * sizeof(real));
   61   tm = (real *) malloc(maxs * sizeof(real));
   62   t0m = array_glb->t0m;
   63   xsm = array_glb->xsm;
   64   zsm = array_glb->zsm;
   65   xrm = array_glb->xrm;
   66   zrm = array_glb->zrm;
   67 
   68   t0m_dim1 = maxr;
   69   zrm_dim1 = maxr;
   70   xrm_dim1 = maxr;
   71   zsm_dim1 = maxr;
   72   xsm_dim1 = maxr;
   73 
   74   count_dim1 = nxmax;
   75   image_dim1 = nxmax;
   76 
   77   count = array_glb->cie[p]->count;
   78   image = array_glb->cie[p]->image;
   79   for (iz = 0; iz < nzmax; ++iz) {
   80     for (ix = 0; ix < nxmax; ++ix) {
   81       image[ix + iz * image_dim1] = (real)0.;
   82       count[ix + iz * count_dim1] = (real)0.;
   83     }
   84   }
   85 /*
   86   if ((p % get_num_proc())  != get_this_proc()) {
   87     printf("%d : Warning %d %d %d\n", get_this_proc(),
   88       p, get_this_proc(), get_num_proc());
   89   }
   90 */
   91 
   92   pi = atan((float)1.) * (float)4.;
   93   pi180 = pi / (float)180.;
   94   /*  master loop runs inversion over common-receiver gathers */
   95   cur = p;
   96   cur_start = 0;
   97   i_2 = rd_glb->ic2;
   98   i_3 = rd_glb->idc;
   99   if (i_3 < 0) {
  100   for (ig = rd_glb->ic1; ig >= i_2; ig += i_3) {
  101     xrt = (float)0.;
  102     xst = (float)0.;
  103     xwavg = (float)0.;
  104     xwmin = (float)1e6;
  105     xwmax = (float)-1e6;
  106     zwmin = (float)1e6;
  107     zwmax = (float)-1e6;
  108     /*   Setup for common-receiver-gathers */
  109     if (direct == 1) {
  110       ir = ig - 1;
  111       ns = 0;
  112       for (is = 0; is < maxs; ++is) {
  113    if (t0m[ir + is * t0m_dim1] == (float)0.) {
  114      goto L2;
  115    }
  116    zw[ns] = zsm[ir + is * zsm_dim1];
  117    zwmin = dmin(zwmin,zw[ns]);
  118    zwmax = dmax(zwmax,zw[ns]);
  119    xw[ns] = xsm[ir + is * xsm_dim1];
  120    xwavg += xw[ns];
  121    xwmin = dmin(xwmin,xw[ns]);
  122    xwmax = dmax(xwmax,xw[ns]);
  123    xrt += xrm[ir + is * xrm_dim1];
  124    tm[ns] = t0m[ir + is * t0m_dim1];
  125    ++ns;
  126 L2:
  127    ;
  128       }
  129       /*    if less than 2 traveltimes, skip to next CRG */
  130       if (ns <= 2) {
  131    goto L500;
  132       }
  133       xwavg /= (real) ns;
  134       xrt /= (real) ns;
  135       zrt = zrm[ir + zrm_dim1];
  136       xg = xrt;
  137       zg = zrt;
  138 #ifdef LOG_STRING
  139       if (wr_glb->frpt == 1) {
  140    fprintf(fp8, "Run#%2d Iter#%2d of%2d BP#%2d of%2d  Rcvr# %3d @(x,z)=%7.1f%7.1f #picks=%3d\n",
  141 
  142        rd_glb->irun, ij, rd_glb->iter, ip, rd_glb->iback, ir+1, xrt, zrt, ns);
  143       }
  144 #endif
  145 #ifdef PRINT_RUN
  146       if (p == 0 /* (rd_glb->gran - 1) */) {
  147         printf("%d %d : Run#%2d Iter#%2d of%2d BP#%2d of%2d  Rcvr# %3d @(x,z)=%7.1f%7.1f #picks=%3d\n",
  148             get_this_proc(), p, rd_glb->irun, ij, rd_glb->iter, ip, rd_glb->iback, ir+1, xrt, zrt, ns);
  149       }
  150 #endif
  151       /*    start/stop angles must be inside control aperture: zwmin to zwmax */
  152       /*    thmax = angle to the lower edge of control aperture */
  153       thmax = atan((zwmax - zrt) / (xwavg - xrt)) / pi180;
  154       /*<   thmin = atan((zwmin - zrt)/(xwavg - xrt))/pi180 >*/
  155       thmin = atan((zwmin - zrt) / (xwavg - xrt)) / pi180;
  156       if (thmin > rd_glb->th2) {
  157    goto L500;
  158       }
  159       if (thmax < rd_glb->th1) {
  160    goto L500;
  161       }
  162       thmin = dmax(thmin - (float)5.,rd_glb->th1);
  163       thmax = dmin(thmax + (float)5.,rd_glb->th2);
  164     }
  165     /*   Setup for common-source-gathers */
  166     if (direct == -1) {
  167       is = ig-1;
  168       ns = 0;
  169       for (ir = 0; ir < maxr; ++ir) {
  170         if (t0m[ir + is * t0m_dim1] == (float)0.) {
  171      goto L3;
  172    }
  173    zw[ns] = zrm[ir + is * zrm_dim1];
  174    zwmin = dmin(zwmin,zw[ns]);
  175    zwmax = dmax(zwmax,zw[ns]);
  176    xw[ns] = xrm[ir + is * xrm_dim1];
  177    xwavg += xw[ns];
  178    xwmin = dmin(xwmin,xw[ns]);
  179    xwmax = dmax(xwmax,xw[ns]);
  180    xst += xsm[ir + is * xsm_dim1];
  181    tm[ns] = t0m[ir + is * t0m_dim1];
  182    ++ns;
  183 L3:
  184    ;
  185       }
  186       /*    if less than 2 traveltimes, skip to next CRG */
  187       if (ns <= 2) {
  188    goto L500;
  189       }
  190       xwavg /= (real) ns;
  191       xst /= (real) ns;
  192       zst = zsm[is * zsm_dim1 + 1];
  193       xg = xst;
  194       zg = zst;
  195 #ifdef LOG_STRING
  196       if (wr_glb->frpt == 1) {
  197    fprintf(fp8, "Run#%2d Iter#%2d of%2d BP#%2d of%2d Source# %3d @(x,z)=%7.1f%7.1f #picks=%3d\n",
  198 
  199        rd_glb->irun, ij, rd_glb->iter, ip, rd_glb->iback, is+1, xst, zst, ns);
  200       }
  201 #endif
  202 #ifdef PRINT_RUN
  203       if (p == 0 /* (rd_glb->gran - 1) */) {
  204         printf("%d %d : Run#%2d Iter#%2d of%2d BP#%2d of%2d Source# %3d @(x,z)=%7.1f%7.1f #picks=%3d\n",
  205             get_this_proc(), p, rd_glb->irun, ij, rd_glb->iter, ip, rd_glb->iback, is+1, xst, zst, ns);
  206       }
  207 #endif
  208       /*    start/stop angles must be inside control aperture: zwmin to zwmax */
  209       /*    thmin = angle to the lower edge of control aperture */
  210       thmin = atan((zwmax - zst) / (xwavg - xst)) / pi180 + (float)180.;
  211       /*    thmax = angle to the upper edge of control aperture */
  212       thmax = atan((zwmin - zst) / (xwavg - xst)) / pi180 + (float)180.;
  213       if (thmin > rd_glb->th2 + (float)180.) {
  214    goto L500;
  215       }
  216       if (thmax < rd_glb->th1 + (float)180.) {
  217    goto L500;
  218       }
  219       thmin = dmax(thmin - (float)5.,rd_glb->th1 + (float)180.);
  220       thmax = dmin(thmax + (float)5.,rd_glb->th2 + (float)180.);
  221     }
  222 #ifdef LOG_STRING
  223     if (wr_glb->frpt == 1) {
  224       fprintf(fp8, "%3d rays in field aperture between  %8.1f %8.1f and %8.1f %8.1f\n",
  225           ns,xwmin,zwmin,xwmax,zwmax);
  226     }
  227 #endif
  228     nray = (r_1 = thmax - thmin, dabs(r_1)) / dabs(rd_glb->thd) + 1;
  229 #ifdef LOG_STRING
  230     if (wr_glb->frpt == 1) {
  231       fprintf(fp8, "%3d synthetic rays launched between angles%7.1f and %7.1f degs\n", nray, thmin, thmax);
  232     }
  233 #endif
  234 
  235 
  236     /*    ***********************************************************
  237 *************/
  238 
  239     icol = 0;
  240     zsimax = (float)-1e6;
  241     zsimin = (float)1e6;
  242     /*   t1 = tan(thmin*pi180) */
  243     /*   t2 = tan(thmax*pi180) */
  244     /*   dt = (t2 - t1)/float(nray-1) */
  245     cur_stop = cur_start + (2 * nray);
  246     while (cur < cur_stop) {
  247         kk = ((cur - cur_start) / nray) + 1;
  248         k = ((cur - cur_start) % nray) + 1;
  249         trad = (thmin + (real) (k - 1) * rd_glb->thd) * pi180;
  250         if (kk == 2) {
  251      trad = -(doublereal)trad;
  252    }
  253    if (kk == 2 && trad == (float)0.) {
  254      goto L400;
  255    }
  256    ncmax = (r_1 = (xwmax - xg) / (rd_glb->dl * cos(trad)), dabs(r_1)) * (float)2.;
  257    rayspace_(s0b, rd_glb->nvb, rd_glb->nhb,
  258        rd_glb->bx0b, rd_glb->bz0b, rd_glb->x0b,
  259        rd_glb->z0b, rd_glb->x1b, rd_glb->z1b, xw, zw,
  260        ns, xwmin, xwmax, ray1, xray1, zray1,
  261        ncmax, xg, zg, trad, rd_glb->dl, &nsamp, &xsi,
  262        &zsi, &rli, rd_glb->isp, direct, &istat);
  263 /*
  264         printf("%d : nsamp %d %lf %lf %lf %d\n", get_this_proc(),
  265           nsamp, xsi, zsi, rli, istat);
  266 */
  267    /*       write out error messages and handle recovery */
  268    if (istat == 1) {
  269      printf("fan # %d ray # %d\n", ig, k);
  270      printf("fatal error: well geometry is not properly defined\n");
  271 #ifdef LOG_STRING
  272      if (wr_glb->frpt == 1) {
  273        fprintf(fp8, "fan # %d ray # %d\n", ig, k);
  274      }
  275      if (wr_glb->frpt == 1) {
  276        fprintf(fp8, "fatal error: well geometry is not properly defined\n");
  277      }
  278 #endif
  279      exit(-1);
  280    }
  281    if (istat == 2) {
  282      /* these lines commented out because they make the reports too large */
  283      /*        write(8,*)'fan #',ig,' ray #',k */
  284      /*        write(8,*)'skip:  ray fell outside defined well geometry' */
  285      goto L400;
  286    }
  287    if (istat == 3) {
  288      /* these lines commented out because they make the reports too large */
  289      /*        write(*,*)'fan #',ig,' ray #',k */
  290      /*        write(*,*)'error: number of raytrace steps > ',ncmax */
  291      /*        if(frpt.eq.1)then */
  292      /*      write(8,*)'fan #',ig,' ray #',k */
  293      /*           write(8,*)'error: */
  294      /*    *      number of raytrace steps >',ncmax */
  295      /*      write(8,*)'Launch angle =',trad/pi180 */
  296      /*      do iw = 1,ns */
  297      /*      write(8,*)iw,xw(iw),zw(iw) */
  298      /*      enddo */
  299      /*        endif */
  300      goto L400;
  301    }
  302    if (istat == 4) {
  303      goto L400;
  304    }
  305    /*    do not use this ray if it falls outside measured aperture */
  306    if (zsi < zwmin || zsi > zwmax) {
  307      goto L400;
  308    }
  309    /*    interpolate raypath */
  310    rayint_(nsamp, xray1, zray1, ray1, rd_glb->idl, xray, zray, rayl);
  311    nstep = (nsamp - 1) * rd_glb->idl + 1;
  312    dl1 = rd_glb->dl / rd_glb->idl;
  313    /*    compute traveltime */
  314    if (ip == 1) {
  315      raytime_(s0, rd_glb->nv, rd_glb->nh,
  316          rd_glb->bx0, rd_glb->bz0, rd_glb->x0,
  317          rd_glb->z0, rd_glb->x1, rd_glb->z1, rd_glb->isp,
  318          nstep, dl1, rayl, xray, zray, time);
  319      tci = time[nstep - 1] * (float)1e3;
  320    } else {
  321      raytime_(s1, nxmax, nzmax, rd_glb->binx,
  322          rd_glb->binz, rd_glb->xmin, rd_glb->zmin,
  323          rd_glb->xmax, rd_glb->zmax, rd_glb->isp,
  324          nstep, dl1, rayl, xray, zray, time);
  325      tci = time[nstep - 1] * (float)1e3;
  326    }
  327    /*    compute the interpolated measured time at the ray/well joint */
  328    tmi = xnterp_(rd_glb->gap, ns, zw, tm, zsi, &istat) * (float)1e3;
  329    if (istat != 0) {
  330      goto L400;
  331    }
  332    /*    do not use this ray if residual is too large */
  333    /*    Residual traveltime */
  334    dta = (r_1 = tmi - tci, dabs(r_1));
  335    if (dta / tmi > bnd2 / (float)100.) {
  336      goto L400;
  337    }
  338    /*          write out ray steps used in inversion */
  339 #ifdef LOG_STRING
  340    if (ig == wr_glb->fray) {
  341      wr_glb->fwrite = 1;
  342      fprintf(fp26, "%d\n", nsamp);
  343      for (ii = 1; ii <= nsamp; ++ii) {
  344        fprintf(fp26, "%f %f\n", xray1[ii - 1], -(doublereal)zray1[ii - 1]);
  345      }
  346    }
  347 #endif
  348    ++icol;
  349    ++ntot;
  350    zsimin = dmin(zsi,zsimin);
  351    zsimax = dmax(zsi,zsimax);
  352    /*    Backproject residuals into image buffers */
  353    resid = (tmi - itom * tci) / rli * cos(trad);
  354    /*    resid = (tmi - itom*tci)/rli */
  355    for (ii = 1; ii <= nstep; ++ii) {
  356      ix = (integer) ((xray[ii - 1] - rd_glb->xmin) / rd_glb->binx);
  357      iz = (integer) ((zray[ii - 1] - rd_glb->zmin) / rd_glb->binz);
  358      if (ix < 0 || ix > (nxmax-1)) {
  359        goto L550;
  360      }
  361      if (iz < 0 || iz > (nzmax-1)) {
  362        goto L550;
  363      }
  364      ++count[ix + iz * count_dim1];
  365      accel = cos(pi * (ii - nstep / 2) / (nstep * (float) 2.)) * (float)1.5;
  366      /* accel = 1.0 */
  367      image[ix + iz * image_dim1] += resid * accel;
  368 L550:
  369      ;
  370    }
  371    /*    calculate traveltime error */
  372    error += dta;
  373    /*    save residuals in microseconds ... for easy display */
  374 
  375 #if defined(TMOD_STRING)
  376    if (wr_glb->fmod == 1 && direct == 1) {
  377      is1 = (zsi - wr_glb->smin) / wr_glb->ds + 1;
  378      if (is1 < 1 || is1 > maxs) {
  379        goto L400;
  380      }
  381      /*  tmod(ig,is1) = (tmi - tci) */
  382      tmod[ig + is1 * tmod_dim1] = tci;
  383    }
  384    if (wr_glb->fmod == 1 && direct == -1) {
  385      ir1 = (zsi - wr_glb->rmin) / wr_glb->dr + 1;
  386      if (ir1 < 1 || ir1 > maxr) {
  387        goto L400;
  388      }
  389      /*  tmod(ir1,ig) = (tmi - tci) */
  390      tmod[ir1 + ig * tmod_dim1] = tci;
  391    }
  392 #endif
  393 L400:
  394    ;
  395         cur += rd_glb->gran;
  396     }
  397     cur_start = cur_stop;
  398 #ifdef LOG_STRING
  399     if (wr_glb->frpt == 1) {
  400       fprintf(fp8, "%3d synthetic rays captured in synth aperture between %8.1f and %8.1f\n\n", icol,zsimin,
  401           zsimax);
  402     }
  403 #endif
  404 L500:
  405     ;
  406   }
  407   } else {
  408    #pragma omp parallel for private (xwmin, ix, iz, thmin, is, ir, xg, ii, ncmax, tci, nstep, zrt, thmax, zsimin, xwmax, xrt, zwmin, ns, cur_stop, zg, zwmax, accel, zsimax, tmi, resid, icol, zst, trad, xwavg, nray, xst, kk, r_1, k, dl1, dta) schedule(static)
  409   for (ig = rd_glb->ic1; ig <= i_2; ig += i_3) {
  410     xrt = (float)0.;
  411     xst = (float)0.;
  412     xwavg = (float)0.;
  413     xwmin = (float)1e6;
  414     xwmax = (float)-1e6;
  415     zwmin = (float)1e6;
  416     zwmax = (float)-1e6;
  417     /*   Setup for common-receiver-gathers */
  418     if (direct == 1) {
  419       ir = ig - 1;
  420       ns = 0;
  421       for (is = 0; is < maxs; ++is) {
  422    if (t0m[ir + is * t0m_dim1] == (float)0.) {
  423      goto L2_;
  424    }
  425    zw[ns] = zsm[ir + is * zsm_dim1];
  426    zwmin = dmin(zwmin,zw[ns]);
  427    zwmax = dmax(zwmax,zw[ns]);
  428    xw[ns] = xsm[ir + is * xsm_dim1];
  429    xwavg += xw[ns];
  430    xwmin = dmin(xwmin,xw[ns]);
  431    xwmax = dmax(xwmax,xw[ns]);
  432    xrt += xrm[ir + is * xrm_dim1];
  433    tm[ns] = t0m[ir + is * t0m_dim1];
  434    ++ns;
  435 L2_:
  436    ;
  437       }
  438       /*    if less than 2 traveltimes, skip to next CRG */
  439       if (ns <= 2) {
  440    goto L500_;
  441       }
  442       xwavg /= (real) ns;
  443       xrt /= (real) ns;
  444       zrt = zrm[ir + zrm_dim1];
  445       xg = xrt;
  446       zg = zrt;
  447 #ifdef LOG_STRING
  448       if (wr_glb->frpt == 1) {
  449    fprintf(fp8, "Run#%2d Iter#%2d of%2d BP#%2d of%2d  Rcvr# %3d @(x,z)=%7.1f%7.1f #picks=%3d\n",
  450 
  451        rd_glb->irun, ij, rd_glb->iter, ip, rd_glb->iback, ir+1, xrt, zrt, ns);
  452       }
  453 #endif
  454 #ifdef PRINT_RUN
  455       if (p == 0 /* (rd_glb->gran - 1) */) {
  456         printf("%d %d : Run#%2d Iter#%2d of%2d BP#%2d of%2d  Rcvr# %3d @(x,z)=%7.1f%7.1f #picks=%3d\n",
  457             get_this_proc(), p, rd_glb->irun, ij, rd_glb->iter, ip, rd_glb->iback, ir+1, xrt, zrt, ns);
  458       }
  459 #endif
  460       /*    start/stop angles must be inside control aperture: zwmin to zwmax */
  461       /*    thmax = angle to the lower edge of control aperture */
  462       thmax = atan((zwmax - zrt) / (xwavg - xrt)) / pi180;
  463       /*<   thmin = atan((zwmin - zrt)/(xwavg - xrt))/pi180 >*/
  464       thmin = atan((zwmin - zrt) / (xwavg - xrt)) / pi180;
  465       if (thmin > rd_glb->th2) {
  466    goto L500_;
  467       }
  468       if (thmax < rd_glb->th1) {
  469    goto L500_;
  470       }
  471       thmin = dmax(thmin - (float)5.,rd_glb->th1);
  472       thmax = dmin(thmax + (float)5.,rd_glb->th2);
  473     }
  474     /*   Setup for common-source-gathers */
  475     if (direct == -1) {
  476       is = ig-1;
  477       ns = 0;
  478       for (ir = 0; ir < maxr; ++ir) {
  479         if (t0m[ir + is * t0m_dim1] == (float)0.) {
  480      goto L3_;
  481    }
  482    zw[ns] = zrm[ir + is * zrm_dim1];
  483    zwmin = dmin(zwmin,zw[ns]);
  484    zwmax = dmax(zwmax,zw[ns]);
  485    xw[ns] = xrm[ir + is * xrm_dim1];
  486    xwavg += xw[ns];
  487    xwmin = dmin(xwmin,xw[ns]);
  488    xwmax = dmax(xwmax,xw[ns]);
  489    xst += xsm[ir + is * xsm_dim1];
  490    tm[ns] = t0m[ir + is * t0m_dim1];
  491    ++ns;
  492 L3_:
  493    ;
  494       }
  495       /*    if less than 2 traveltimes, skip to next CRG */
  496       if (ns <= 2) {
  497    goto L500_;
  498       }
  499       xwavg /= (real) ns;
  500       xst /= (real) ns;
  501       zst = zsm[is * zsm_dim1 + 1];
  502       xg = xst;
  503       zg = zst;
  504 #ifdef LOG_STRING
  505       if (wr_glb->frpt == 1) {
  506    fprintf(fp8, "Run#%2d Iter#%2d of%2d BP#%2d of%2d Source# %3d @(x,z)=%7.1f%7.1f #picks=%3d\n",
  507 
  508        rd_glb->irun, ij, rd_glb->iter, ip, rd_glb->iback, is+1, xst, zst, ns);
  509       }
  510 #endif
  511 #ifdef PRINT_RUN
  512       if (p == 0 /* (rd_glb->gran - 1) */) {
  513         printf("%d %d : Run#%2d Iter#%2d of%2d BP#%2d of%2d Source# %3d @(x,z)=%7.1f%7.1f #picks=%3d\n",
  514             get_this_proc(), p, rd_glb->irun, ij, rd_glb->iter, ip, rd_glb->iback, is+1, xst, zst, ns);
  515       }
  516 #endif
  517       /*    start/stop angles must be inside control aperture: zwmin to zwmax */
  518       /*    thmin = angle to the lower edge of control aperture */
  519       thmin = atan((zwmax - zst) / (xwavg - xst)) / pi180 + (float)180.;
  520       /*    thmax = angle to the upper edge of control aperture */
  521       thmax = atan((zwmin - zst) / (xwavg - xst)) / pi180 + (float)180.;
  522       if (thmin > rd_glb->th2 + (float)180.) {
  523    goto L500_;
  524       }
  525       if (thmax < rd_glb->th1 + (float)180.) {
  526    goto L500_;
  527       }
  528       thmin = dmax(thmin - (float)5.,rd_glb->th1 + (float)180.);
  529       thmax = dmin(thmax + (float)5.,rd_glb->th2 + (float)180.);
  530     }
  531 #ifdef LOG_STRING
  532     if (wr_glb->frpt == 1) {
  533       fprintf(fp8, "%3d rays in field aperture between  %8.1f %8.1f and %8.1f %8.1f\n",
  534           ns,xwmin,zwmin,xwmax,zwmax);
  535     }
  536 #endif
  537     nray = (r_1 = thmax - thmin, dabs(r_1)) / dabs(rd_glb->thd) + 1;
  538 #ifdef LOG_STRING
  539     if (wr_glb->frpt == 1) {
  540       fprintf(fp8, "%3d synthetic rays launched between angles%7.1f and %7.1f degs\n", nray, thmin, thmax);
  541     }
  542 #endif
  543 
  544 
  545     /*    ***********************************************************
  546 *************/
  547 
  548     icol = 0;
  549     zsimax = (float)-1e6;
  550     zsimin = (float)1e6;
  551     /*   t1 = tan(thmin*pi180) */
  552     /*   t2 = tan(thmax*pi180) */
  553     /*   dt = (t2 - t1)/float(nray-1) */
  554     cur_stop = cur_start + (2 * nray);
  555     while (cur < cur_stop) {
  556         kk = ((cur - cur_start) / nray) + 1;
  557         k = ((cur - cur_start) % nray) + 1;
  558         trad = (thmin + (real) (k - 1) * rd_glb->thd) * pi180;
  559         if (kk == 2) {
  560      trad = -(doublereal)trad;
  561    }
  562    if (kk == 2 && trad == (float)0.) {
  563      goto L400_;
  564    }
  565    ncmax = (r_1 = (xwmax - xg) / (rd_glb->dl * cos(trad)), dabs(r_1)) * (float)2.;
  566    rayspace_(s0b, rd_glb->nvb, rd_glb->nhb,
  567        rd_glb->bx0b, rd_glb->bz0b, rd_glb->x0b,
  568        rd_glb->z0b, rd_glb->x1b, rd_glb->z1b, xw, zw,
  569        ns, xwmin, xwmax, ray1, xray1, zray1,
  570        ncmax, xg, zg, trad, rd_glb->dl, &nsamp, &xsi,
  571        &zsi, &rli, rd_glb->isp, direct, &istat);
  572 /*
  573         printf("%d : nsamp %d %lf %lf %lf %d\n", get_this_proc(),
  574           nsamp, xsi, zsi, rli, istat);
  575 */
  576    /*       write out error messages and handle recovery */
  577    if (istat == 1) {
  578      printf("fan # %d ray # %d\n", ig, k);
  579      printf("fatal error: well geometry is not properly defined\n");
  580 #ifdef LOG_STRING
  581      if (wr_glb->frpt == 1) {
  582        fprintf(fp8, "fan # %d ray # %d\n", ig, k);
  583      }
  584      if (wr_glb->frpt == 1) {
  585        fprintf(fp8, "fatal error: well geometry is not properly defined\n");
  586      }
  587 #endif
  588      exit(-1);
  589    }
  590    if (istat == 2) {
  591      /* these lines commented out because they make the reports too large */
  592      /*        write(8,*)'fan #',ig,' ray #',k */
  593      /*        write(8,*)'skip:  ray fell outside defined well geometry' */
  594      goto L400_;
  595    }
  596    if (istat == 3) {
  597      /* these lines commented out because they make the reports too large */
  598      /*        write(*,*)'fan #',ig,' ray #',k */
  599      /*        write(*,*)'error: number of raytrace steps > ',ncmax */
  600      /*        if(frpt.eq.1)then */
  601      /*      write(8,*)'fan #',ig,' ray #',k */
  602      /*           write(8,*)'error: */
  603      /*    *      number of raytrace steps >',ncmax */
  604      /*      write(8,*)'Launch angle =',trad/pi180 */
  605      /*      do iw = 1,ns */
  606      /*      write(8,*)iw,xw(iw),zw(iw) */
  607      /*      enddo */
  608      /*        endif */
  609      goto L400_;
  610    }
  611    if (istat == 4) {
  612      goto L400_;
  613    }
  614    /*    do not use this ray if it falls outside measured aperture */
  615    if (zsi < zwmin || zsi > zwmax) {
  616      goto L400_;
  617    }
  618    /*    interpolate raypath */
  619    rayint_(nsamp, xray1, zray1, ray1, rd_glb->idl, xray, zray, rayl);
  620    nstep = (nsamp - 1) * rd_glb->idl + 1;
  621    dl1 = rd_glb->dl / rd_glb->idl;
  622    /*    compute traveltime */
  623    if (ip == 1) {
  624      raytime_(s0, rd_glb->nv, rd_glb->nh,
  625          rd_glb->bx0, rd_glb->bz0, rd_glb->x0,
  626          rd_glb->z0, rd_glb->x1, rd_glb->z1, rd_glb->isp,
  627          nstep, dl1, rayl, xray, zray, time);
  628      tci = time[nstep - 1] * (float)1e3;
  629    } else {
  630      raytime_(s1, nxmax, nzmax, rd_glb->binx,
  631          rd_glb->binz, rd_glb->xmin, rd_glb->zmin,
  632          rd_glb->xmax, rd_glb->zmax, rd_glb->isp,
  633          nstep, dl1, rayl, xray, zray, time);
  634      tci = time[nstep - 1] * (float)1e3;
  635    }
  636    /*    compute the interpolated measured time at the ray/well joint */
  637    tmi = xnterp_(rd_glb->gap, ns, zw, tm, zsi, &istat) * (float)1e3;
  638    if (istat != 0) {
  639      goto L400_;
  640    }
  641    /*    do not use this ray if residual is too large */
  642    /*    Residual traveltime */
  643    dta = (r_1 = tmi - tci, dabs(r_1));
  644    if (dta / tmi > bnd2 / (float)100.) {
  645      goto L400_;
  646    }
  647    /*          write out ray steps used in inversion */
  648 #ifdef LOG_STRING
  649    if (ig == wr_glb->fray) {
  650      wr_glb->fwrite = 1;
  651      fprintf(fp26, "%d\n", nsamp);
  652      for (ii = 1; ii <= nsamp; ++ii) {
  653        fprintf(fp26, "%f %f\n", xray1[ii - 1], -(doublereal)zray1[ii - 1]);
  654      }
  655    }
  656 #endif
  657    ++icol;
  658    ++ntot;
  659    zsimin = dmin(zsi,zsimin);
  660    zsimax = dmax(zsi,zsimax);
  661    /*    Backproject residuals into image buffers */
  662    resid = (tmi - itom * tci) / rli * cos(trad);
  663    /*    resid = (tmi - itom*tci)/rli */
  664    for (ii = 1; ii <= nstep; ++ii) {
  665      ix = (integer) ((xray[ii - 1] - rd_glb->xmin) / rd_glb->binx);
  666      iz = (integer) ((zray[ii - 1] - rd_glb->zmin) / rd_glb->binz);
  667      if (ix < 0 || ix > (nxmax-1)) {
  668        goto L550_;
  669      }
  670      if (iz < 0 || iz > (nzmax-1)) {
  671        goto L550_;
  672      }
  673      ++count[ix + iz * count_dim1];
  674      accel = cos(pi * (ii - nstep / 2) / (nstep * (float) 2.)) * (float)1.5;
  675      /* accel = 1.0 */
  676      image[ix + iz * image_dim1] += resid * accel;
  677 L550_:
  678      ;
  679    }
  680    /*    calculate traveltime error */
  681    error += dta;
  682    /*    save residuals in microseconds ... for easy display */
  683 
  684 #if defined(TMOD_STRING)
  685    if (wr_glb->fmod == 1 && direct == 1) {
  686      is1 = (zsi - wr_glb->smin) / wr_glb->ds + 1;
  687      if (is1 < 1 || is1 > maxs) {
  688        goto L400_;
  689      }
  690      /*  tmod(ig,is1) = (tmi - tci) */
  691      tmod[ig + is1 * tmod_dim1] = tci;
  692    }
  693    if (wr_glb->fmod == 1 && direct == -1) {
  694      ir1 = (zsi - wr_glb->rmin) / wr_glb->dr + 1;
  695      if (ir1 < 1 || ir1 > maxr) {
  696        goto L400_;
  697      }
  698      /*  tmod(ir1,ig) = (tmi - tci) */
  699      tmod[ir1 + ig * tmod_dim1] = tci;
  700    }
  701 #endif
  702 L400_:
  703    ;
  704         cur += rd_glb->gran;
  705     }
  706     cur_start = cur_stop;
  707 #ifdef LOG_STRING
  708     if (wr_glb->frpt == 1) {
  709       fprintf(fp8, "%3d synthetic rays captured in synth aperture between %8.1f and %8.1f\n\n", icol,zsimin,
  710           zsimax);
  711     }
  712 #endif
  713 L500_:
  714     ;
  715   }
  716   }
  717   array_glb->cie[p]->error = error;
  718   free(xw);
  719   free(zw);
  720   free(tm);
  721   free(xray1);
  722   for (i = 0, rp = (p | (1 << i)); (rp < rd_glb->gran) && (rp != p); i++, rp = (p | (1 << i))) {
  723     array_glb->cie[p]->error += array_glb->cie[rp]->error;
  724     rc = array_glb->cie[rp]->count;
  725     ri = array_glb->cie[rp]->image;
  726     for (iz = 0; iz < nzmax; ++iz) {
  727       for (ix = 0; ix < nxmax; ++ix) {
  728         image[ix + iz * image_dim1] += ri[ix + iz * image_dim1];
  729         count[ix + iz * count_dim1] += rc[ix + iz * count_dim1];
  730       }
  731     }
  732   }
  733 }
  734 
  735 doloop(ij, ip, bnd2, itom, nxmax, nzmax, maxs, maxr, direct,
  736   sh_s0, sh_s0b, sh_s1, sh_rd_glb, sh_array_glb)
  737 integer ij, ip, itom, nxmax, nzmax, maxs, maxr;
  738 integer direct;
  739 double bnd2;
  740 sh_ptr_real sh_s0, sh_s0b, sh_s1;
  741 sh_ptr_rd_glb sh_rd_glb;
  742 sh_ptr_array_glb sh_array_glb;
  743 {
  744   integer count_dim1, image_dim1, p, rp, i, ix, iz;
  745   rd_wr_ptr_real count, image;
  746   rd_ptr_rd_glb rd_glb = sh_rd_glb;
  747   rd_ptr_array_glb array_glb = sh_array_glb;
  748 
  749   count_dim1 = nxmax;
  750   image_dim1 = nxmax;
  751 
  752   for (p = rd_glb->gran - 1; p >= 0; p--) {
  753       runloop(p, ij, ip, bnd2, itom, nxmax, nzmax, maxs, maxr,
  754               direct, sh_s0, sh_s0b, sh_s1, sh_rd_glb, sh_array_glb);
  755   }
  756 }