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 }