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   for (ig = rd_glb->ic1; ig <= i_2; ig += i_3) {
  409     xrt = (float)0.;
  410     xst = (float)0.;
  411     xwavg = (float)0.;
  412     xwmin = (float)1e6;
  413     xwmax = (float)-1e6;
  414     zwmin = (float)1e6;
  415     zwmax = (float)-1e6;
  416     /*   Setup for common-receiver-gathers */
  417     if (direct == 1) {
  418       ir = ig - 1;
  419       ns = 0;
  420       for (is = 0; is < maxs; ++is) {
  421    if (t0m[ir + is * t0m_dim1] == (float)0.) {
  422      goto L2_;
  423    }
  424    zw[ns] = zsm[ir + is * zsm_dim1];
  425    zwmin = dmin(zwmin,zw[ns]);
  426    zwmax = dmax(zwmax,zw[ns]);
  427    xw[ns] = xsm[ir + is * xsm_dim1];
  428    xwavg += xw[ns];
  429    xwmin = dmin(xwmin,xw[ns]);
  430    xwmax = dmax(xwmax,xw[ns]);
  431    xrt += xrm[ir + is * xrm_dim1];
  432    tm[ns] = t0m[ir + is * t0m_dim1];
  433    ++ns;
  434 L2_:
  435    ;
  436       }
  437       /*    if less than 2 traveltimes, skip to next CRG */
  438       if (ns <= 2) {
  439    goto L500_;
  440       }
  441       xwavg /= (real) ns;
  442       xrt /= (real) ns;
  443       zrt = zrm[ir + zrm_dim1];
  444       xg = xrt;
  445       zg = zrt;
  446 #ifdef LOG_STRING
  447       if (wr_glb->frpt == 1) {
  448    fprintf(fp8, "Run#%2d Iter#%2d of%2d BP#%2d of%2d  Rcvr# %3d @(x,z)=%7.1f%7.1f #picks=%3d\n",
  449 
  450        rd_glb->irun, ij, rd_glb->iter, ip, rd_glb->iback, ir+1, xrt, zrt, ns);
  451       }
  452 #endif
  453 #ifdef PRINT_RUN
  454       if (p == 0 /* (rd_glb->gran - 1) */) {
  455         printf("%d %d : Run#%2d Iter#%2d of%2d BP#%2d of%2d  Rcvr# %3d @(x,z)=%7.1f%7.1f #picks=%3d\n",
  456             get_this_proc(), p, rd_glb->irun, ij, rd_glb->iter, ip, rd_glb->iback, ir+1, xrt, zrt, ns);
  457       }
  458 #endif
  459       /*    start/stop angles must be inside control aperture: zwmin to zwmax */
  460       /*    thmax = angle to the lower edge of control aperture */
  461       thmax = atan((zwmax - zrt) / (xwavg - xrt)) / pi180;
  462       /*<   thmin = atan((zwmin - zrt)/(xwavg - xrt))/pi180 >*/
  463       thmin = atan((zwmin - zrt) / (xwavg - xrt)) / pi180;
  464       if (thmin > rd_glb->th2) {
  465    goto L500_;
  466       }
  467       if (thmax < rd_glb->th1) {
  468    goto L500_;
  469       }
  470       thmin = dmax(thmin - (float)5.,rd_glb->th1);
  471       thmax = dmin(thmax + (float)5.,rd_glb->th2);
  472     }
  473     /*   Setup for common-source-gathers */
  474     if (direct == -1) {
  475       is = ig-1;
  476       ns = 0;
  477       for (ir = 0; ir < maxr; ++ir) {
  478         if (t0m[ir + is * t0m_dim1] == (float)0.) {
  479      goto L3_;
  480    }
  481    zw[ns] = zrm[ir + is * zrm_dim1];
  482    zwmin = dmin(zwmin,zw[ns]);
  483    zwmax = dmax(zwmax,zw[ns]);
  484    xw[ns] = xrm[ir + is * xrm_dim1];
  485    xwavg += xw[ns];
  486    xwmin = dmin(xwmin,xw[ns]);
  487    xwmax = dmax(xwmax,xw[ns]);
  488    xst += xsm[ir + is * xsm_dim1];
  489    tm[ns] = t0m[ir + is * t0m_dim1];
  490    ++ns;
  491 L3_:
  492    ;
  493       }
  494       /*    if less than 2 traveltimes, skip to next CRG */
  495       if (ns <= 2) {
  496    goto L500_;
  497       }
  498       xwavg /= (real) ns;
  499       xst /= (real) ns;
  500       zst = zsm[is * zsm_dim1 + 1];
  501       xg = xst;
  502       zg = zst;
  503 #ifdef LOG_STRING
  504       if (wr_glb->frpt == 1) {
  505    fprintf(fp8, "Run#%2d Iter#%2d of%2d BP#%2d of%2d Source# %3d @(x,z)=%7.1f%7.1f #picks=%3d\n",
  506 
  507        rd_glb->irun, ij, rd_glb->iter, ip, rd_glb->iback, is+1, xst, zst, ns);
  508       }
  509 #endif
  510 #ifdef PRINT_RUN
  511       if (p == 0 /* (rd_glb->gran - 1) */) {
  512         printf("%d %d : Run#%2d Iter#%2d of%2d BP#%2d of%2d Source# %3d @(x,z)=%7.1f%7.1f #picks=%3d\n",
  513             get_this_proc(), p, rd_glb->irun, ij, rd_glb->iter, ip, rd_glb->iback, is+1, xst, zst, ns);
  514       }
  515 #endif
  516       /*    start/stop angles must be inside control aperture: zwmin to zwmax */
  517       /*    thmin = angle to the lower edge of control aperture */
  518       thmin = atan((zwmax - zst) / (xwavg - xst)) / pi180 + (float)180.;
  519       /*    thmax = angle to the upper edge of control aperture */
  520       thmax = atan((zwmin - zst) / (xwavg - xst)) / pi180 + (float)180.;
  521       if (thmin > rd_glb->th2 + (float)180.) {
  522    goto L500_;
  523       }
  524       if (thmax < rd_glb->th1 + (float)180.) {
  525    goto L500_;
  526       }
  527       thmin = dmax(thmin - (float)5.,rd_glb->th1 + (float)180.);
  528       thmax = dmin(thmax + (float)5.,rd_glb->th2 + (float)180.);
  529     }
  530 #ifdef LOG_STRING
  531     if (wr_glb->frpt == 1) {
  532       fprintf(fp8, "%3d rays in field aperture between  %8.1f %8.1f and %8.1f %8.1f\n",
  533           ns,xwmin,zwmin,xwmax,zwmax);
  534     }
  535 #endif
  536     nray = (r_1 = thmax - thmin, dabs(r_1)) / dabs(rd_glb->thd) + 1;
  537 #ifdef LOG_STRING
  538     if (wr_glb->frpt == 1) {
  539       fprintf(fp8, "%3d synthetic rays launched between angles%7.1f and %7.1f degs\n", nray, thmin, thmax);
  540     }
  541 #endif
  542 
  543 
  544     /*    ***********************************************************
  545 *************/
  546 
  547     icol = 0;
  548     zsimax = (float)-1e6;
  549     zsimin = (float)1e6;
  550     /*   t1 = tan(thmin*pi180) */
  551     /*   t2 = tan(thmax*pi180) */
  552     /*   dt = (t2 - t1)/float(nray-1) */
  553     cur_stop = cur_start + (2 * nray);
  554     while (cur < cur_stop) {
  555         kk = ((cur - cur_start) / nray) + 1;
  556         k = ((cur - cur_start) % nray) + 1;
  557         trad = (thmin + (real) (k - 1) * rd_glb->thd) * pi180;
  558         if (kk == 2) {
  559      trad = -(doublereal)trad;
  560    }
  561    if (kk == 2 && trad == (float)0.) {
  562      goto L400_;
  563    }
  564    ncmax = (r_1 = (xwmax - xg) / (rd_glb->dl * cos(trad)), dabs(r_1)) * (float)2.;
  565    rayspace_(s0b, rd_glb->nvb, rd_glb->nhb,
  566        rd_glb->bx0b, rd_glb->bz0b, rd_glb->x0b,
  567        rd_glb->z0b, rd_glb->x1b, rd_glb->z1b, xw, zw,
  568        ns, xwmin, xwmax, ray1, xray1, zray1,
  569        ncmax, xg, zg, trad, rd_glb->dl, &nsamp, &xsi,
  570        &zsi, &rli, rd_glb->isp, direct, &istat);
  571 /*
  572         printf("%d : nsamp %d %lf %lf %lf %d\n", get_this_proc(),
  573           nsamp, xsi, zsi, rli, istat);
  574 */
  575    /*       write out error messages and handle recovery */
  576    if (istat == 1) {
  577      printf("fan # %d ray # %d\n", ig, k);
  578      printf("fatal error: well geometry is not properly defined\n");
  579 #ifdef LOG_STRING
  580      if (wr_glb->frpt == 1) {
  581        fprintf(fp8, "fan # %d ray # %d\n", ig, k);
  582      }
  583      if (wr_glb->frpt == 1) {
  584        fprintf(fp8, "fatal error: well geometry is not properly defined\n");
  585      }
  586 #endif
  587      exit(-1);
  588    }
  589    if (istat == 2) {
  590      /* these lines commented out because they make the reports too large */
  591      /*        write(8,*)'fan #',ig,' ray #',k */
  592      /*        write(8,*)'skip:  ray fell outside defined well geometry' */
  593      goto L400_;
  594    }
  595    if (istat == 3) {
  596      /* these lines commented out because they make the reports too large */
  597      /*        write(*,*)'fan #',ig,' ray #',k */
  598      /*        write(*,*)'error: number of raytrace steps > ',ncmax */
  599      /*        if(frpt.eq.1)then */
  600      /*      write(8,*)'fan #',ig,' ray #',k */
  601      /*           write(8,*)'error: */
  602      /*    *      number of raytrace steps >',ncmax */
  603      /*      write(8,*)'Launch angle =',trad/pi180 */
  604      /*      do iw = 1,ns */
  605      /*      write(8,*)iw,xw(iw),zw(iw) */
  606      /*      enddo */
  607      /*        endif */
  608      goto L400_;
  609    }
  610    if (istat == 4) {
  611      goto L400_;
  612    }
  613    /*    do not use this ray if it falls outside measured aperture */
  614    if (zsi < zwmin || zsi > zwmax) {
  615      goto L400_;
  616    }
  617    /*    interpolate raypath */
  618    rayint_(nsamp, xray1, zray1, ray1, rd_glb->idl, xray, zray, rayl);
  619    nstep = (nsamp - 1) * rd_glb->idl + 1;
  620    dl1 = rd_glb->dl / rd_glb->idl;
  621    /*    compute traveltime */
  622    if (ip == 1) {
  623      raytime_(s0, rd_glb->nv, rd_glb->nh,
  624          rd_glb->bx0, rd_glb->bz0, rd_glb->x0,
  625          rd_glb->z0, rd_glb->x1, rd_glb->z1, rd_glb->isp,
  626          nstep, dl1, rayl, xray, zray, time);
  627      tci = time[nstep - 1] * (float)1e3;
  628    } else {
  629      raytime_(s1, nxmax, nzmax, rd_glb->binx,
  630          rd_glb->binz, rd_glb->xmin, rd_glb->zmin,
  631          rd_glb->xmax, rd_glb->zmax, rd_glb->isp,
  632          nstep, dl1, rayl, xray, zray, time);
  633      tci = time[nstep - 1] * (float)1e3;
  634    }
  635    /*    compute the interpolated measured time at the ray/well joint */
  636    tmi = xnterp_(rd_glb->gap, ns, zw, tm, zsi, &istat) * (float)1e3;
  637    if (istat != 0) {
  638      goto L400_;
  639    }
  640    /*    do not use this ray if residual is too large */
  641    /*    Residual traveltime */
  642    dta = (r_1 = tmi - tci, dabs(r_1));
  643    if (dta / tmi > bnd2 / (float)100.) {
  644      goto L400_;
  645    }
  646    /*          write out ray steps used in inversion */
  647 #ifdef LOG_STRING
  648    if (ig == wr_glb->fray) {
  649      wr_glb->fwrite = 1;
  650      fprintf(fp26, "%d\n", nsamp);
  651      for (ii = 1; ii <= nsamp; ++ii) {
  652        fprintf(fp26, "%f %f\n", xray1[ii - 1], -(doublereal)zray1[ii - 1]);
  653      }
  654    }
  655 #endif
  656    ++icol;
  657    ++ntot;
  658    zsimin = dmin(zsi,zsimin);
  659    zsimax = dmax(zsi,zsimax);
  660    /*    Backproject residuals into image buffers */
  661    resid = (tmi - itom * tci) / rli * cos(trad);
  662    /*    resid = (tmi - itom*tci)/rli */
  663    for (ii = 1; ii <= nstep; ++ii) {
  664      ix = (integer) ((xray[ii - 1] - rd_glb->xmin) / rd_glb->binx);
  665      iz = (integer) ((zray[ii - 1] - rd_glb->zmin) / rd_glb->binz);
  666      if (ix < 0 || ix > (nxmax-1)) {
  667        goto L550_;
  668      }
  669      if (iz < 0 || iz > (nzmax-1)) {
  670        goto L550_;
  671      }
  672      ++count[ix + iz * count_dim1];
  673      accel = cos(pi * (ii - nstep / 2) / (nstep * (float) 2.)) * (float)1.5;
  674      /* accel = 1.0 */
  675      image[ix + iz * image_dim1] += resid * accel;
  676 L550_:
  677      ;
  678    }
  679    /*    calculate traveltime error */
  680    error += dta;
  681    /*    save residuals in microseconds ... for easy display */
  682 
  683 #if defined(TMOD_STRING)
  684    if (wr_glb->fmod == 1 && direct == 1) {
  685      is1 = (zsi - wr_glb->smin) / wr_glb->ds + 1;
  686      if (is1 < 1 || is1 > maxs) {
  687        goto L400_;
  688      }
  689      /*  tmod(ig,is1) = (tmi - tci) */
  690      tmod[ig + is1 * tmod_dim1] = tci;
  691    }
  692    if (wr_glb->fmod == 1 && direct == -1) {
  693      ir1 = (zsi - wr_glb->rmin) / wr_glb->dr + 1;
  694      if (ir1 < 1 || ir1 > maxr) {
  695        goto L400_;
  696      }
  697      /*  tmod(ir1,ig) = (tmi - tci) */
  698      tmod[ir1 + ig * tmod_dim1] = tci;
  699    }
  700 #endif
  701 L400_:
  702    ;
  703         cur += rd_glb->gran;
  704     }
  705     cur_start = cur_stop;
  706 #ifdef LOG_STRING
  707     if (wr_glb->frpt == 1) {
  708       fprintf(fp8, "%3d synthetic rays captured in synth aperture between %8.1f and %8.1f\n\n", icol,zsimin,
  709           zsimax);
  710     }
  711 #endif
  712 L500_:
  713     ;
  714   }
  715   }
  716   array_glb->cie[p]->error = error;
  717   free(xw);
  718   free(zw);
  719   free(tm);
  720   free(xray1);
  721   for (i = 0, rp = (p | (1 << i)); (rp < rd_glb->gran) && (rp != p); i++, rp = (p | (1 << i))) {
  722     array_glb->cie[p]->error += array_glb->cie[rp]->error;
  723     rc = array_glb->cie[rp]->count;
  724     ri = array_glb->cie[rp]->image;
  725     for (iz = 0; iz < nzmax; ++iz) {
  726       for (ix = 0; ix < nxmax; ++ix) {
  727         image[ix + iz * image_dim1] += ri[ix + iz * image_dim1];
  728         count[ix + iz * count_dim1] += rc[ix + iz * count_dim1];
  729       }
  730     }
  731   }
  732 }
  733 
  734 doloop(ij, ip, bnd2, itom, nxmax, nzmax, maxs, maxr, direct,
  735   sh_s0, sh_s0b, sh_s1, sh_rd_glb, sh_array_glb)
  736 integer ij, ip, itom, nxmax, nzmax, maxs, maxr;
  737 integer direct;
  738 double bnd2;
  739 sh_ptr_real sh_s0, sh_s0b, sh_s1;
  740 sh_ptr_rd_glb sh_rd_glb;
  741 sh_ptr_array_glb sh_array_glb;
  742 {
  743   integer count_dim1, image_dim1, p, rp, i, ix, iz;
  744   rd_wr_ptr_real count, image;
  745   rd_ptr_rd_glb rd_glb = sh_rd_glb;
  746   rd_ptr_array_glb array_glb = sh_array_glb;
  747 
  748   count_dim1 = nxmax;
  749   image_dim1 = nxmax;
  750 
  751    #pragma omp parallel for schedule(static)
  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 }