|
| 1 | +/* |
| 2 | + *************************************************************************************** |
| 3 | + * |
| 4 | + * NAME: |
| 5 | + * grd3d_get_randomline.c |
| 6 | + * |
| 7 | + * AUTHOR(S): |
| 8 | + * Jan C. Rivenaes |
| 9 | + * |
| 10 | + * DESCRIPTION: |
| 11 | + * Given X Y Z vectors, return a a randomline array from a 3D grid property |
| 12 | + * |
| 13 | + * ARGUMENTS: |
| 14 | + * xvec, yvec i Arrays coords XY |
| 15 | + * zmin, zmax i Vertical range |
| 16 | + * nzsam i Vertical sampling numbering |
| 17 | + * mcol, mrow i Number of rows/cols for maps |
| 18 | + * xori..rotation i Map settings |
| 19 | + * maptopi..mapbasj i Map arrays for I J top/base |
| 20 | + * nx ny nz i Grid dimensions |
| 21 | + * p_zcorn_v i Grid Zcorn |
| 22 | + * p_coord_v i Grid ZCORN |
| 23 | + * p_acnum_v i Grid ACTNUM |
| 24 | + * p_val_v i 3D Grid values |
| 25 | + * p_zcornone_v i Grid ZCORN |
| 26 | + * p_acnumone_v i Grid ACTNUM |
| 27 | + * value o Randomline array |
| 28 | + * option i For later |
| 29 | + * debug i Debug level |
| 30 | + * |
| 31 | + * RETURNS: |
| 32 | + * Array length, -1 if fail |
| 33 | + * |
| 34 | + * TODO/ISSUES/BUGS: |
| 35 | + * |
| 36 | + * LICENCE: |
| 37 | + * cf. XTGeo LICENSE |
| 38 | + *************************************************************************************** |
| 39 | + */ |
| 40 | + |
| 41 | + |
| 42 | +#include "libxtg.h" |
| 43 | +#include "libxtg_.h" |
| 44 | + |
| 45 | + |
| 46 | +/* |
| 47 | +**************************************************************************************** |
| 48 | +* private function |
| 49 | +**************************************************************************************** |
| 50 | +*/ |
| 51 | + |
| 52 | +void _get_ij_range(int *i1, int *i2, int *j1, int *j2, double xc, double yc, int mcol, |
| 53 | + int mrow, double xori, double yori, double xinc, double yinc, |
| 54 | + int yflip, double rotation, double *maptopi, double *maptopj, |
| 55 | + double *mapbasi, double *mapbasj, int debug) |
| 56 | +{ |
| 57 | + char sbn[24] = "_get_ijrange"; |
| 58 | + long nmap; |
| 59 | + int itop, jtop, ibas, jbas, ii1, ii2, jj1, jj2; |
| 60 | + |
| 61 | + xtgverbose(debug); |
| 62 | + |
| 63 | + nmap = mcol * mrow; |
| 64 | + |
| 65 | + /* get map value for I J from x y */ |
| 66 | + itop = surf_get_z_from_xy(xc, yc, mcol, mrow, xori, yori, xinc, yinc, |
| 67 | + yflip, rotation, maptopi, nmap, debug); |
| 68 | + jtop = surf_get_z_from_xy(xc, yc, mcol, mrow, xori, yori, xinc, yinc, |
| 69 | + yflip, rotation, maptopj, nmap, debug); |
| 70 | + ibas = surf_get_z_from_xy(xc, yc, mcol, mrow, xori, yori, xinc, yinc, |
| 71 | + yflip, rotation, mapbasi, nmap, debug); |
| 72 | + jbas = surf_get_z_from_xy(xc, yc, mcol, mrow, xori, yori, xinc, yinc, |
| 73 | + yflip, rotation, mapbasj, nmap, debug); |
| 74 | + |
| 75 | + if (debug > 1) xtg_speak(sbn, 2, "ITOP IBAS JTOP JBAS %d %d %d %d...", |
| 76 | + itop, ibas, jtop, jbas); |
| 77 | + |
| 78 | + if (itop <= ibas){ |
| 79 | + ii1 = itop; |
| 80 | + ii2 = ibas; |
| 81 | + } |
| 82 | + else { |
| 83 | + ii1 = ibas; |
| 84 | + ii2 = itop; |
| 85 | + } |
| 86 | + |
| 87 | + /* extend with one to be sure */ |
| 88 | + if (ii1 > 1) ii1--; |
| 89 | + if (ii2 < mcol) ii2++; |
| 90 | + |
| 91 | + if (jtop <= jbas){ |
| 92 | + jj1 = jtop; |
| 93 | + jj2 = jbas; |
| 94 | + } |
| 95 | + else { |
| 96 | + jj1 = jbas; |
| 97 | + jj2 = jtop; |
| 98 | + } |
| 99 | + |
| 100 | + /* extend with one to be sure */ |
| 101 | + if (jj1 > 1) jj1--; |
| 102 | + if (jj2 < mrow) jj2++; |
| 103 | + |
| 104 | + *i1 = ii1; |
| 105 | + *i2 = ii2; |
| 106 | + *j1 = jj1; |
| 107 | + *j2 = jj2; |
| 108 | + |
| 109 | +} |
| 110 | + |
| 111 | +/* |
| 112 | +**************************************************************************************** |
| 113 | +* public function |
| 114 | +**************************************************************************************** |
| 115 | +*/ |
| 116 | + |
| 117 | +int grd3d_get_randomline( |
| 118 | + double *xvec, |
| 119 | + long nxvec, |
| 120 | + double *yvec, |
| 121 | + long nyvec, |
| 122 | + double zmin, |
| 123 | + double zmax, |
| 124 | + int nzsam, |
| 125 | + |
| 126 | + int mcol, |
| 127 | + int mrow, |
| 128 | + double xori, |
| 129 | + double yori, |
| 130 | + double xinc, |
| 131 | + double yinc, |
| 132 | + double rotation, |
| 133 | + int yflip, |
| 134 | + double *maptopi, |
| 135 | + double *maptopj, |
| 136 | + double *mapbasi, |
| 137 | + double *mapbasj, |
| 138 | + |
| 139 | + int nx, |
| 140 | + int ny, |
| 141 | + int nz, |
| 142 | + double *p_coor_v, |
| 143 | + double *p_zcorn_v, |
| 144 | + int *p_actnum_v, |
| 145 | + double *p_val_v, |
| 146 | + double *p_zcornone_v, |
| 147 | + int *p_actnumone_v, |
| 148 | + |
| 149 | + double *values, |
| 150 | + long nvalues, |
| 151 | + |
| 152 | + int option, |
| 153 | + int debug |
| 154 | + ) |
| 155 | +{ |
| 156 | + /* locals */ |
| 157 | + char sbn[24] = "grd3d_get_randomline"; |
| 158 | + int ib, ic, izc, ier, ios, i1, i2, j1, j2, k1, k2; |
| 159 | + long ibs1, ibs2; |
| 160 | + double zsam, xc, yc, zc; |
| 161 | + double value, *p_dummy_v; |
| 162 | + |
| 163 | + xtgverbose(debug); |
| 164 | + |
| 165 | + if (debug > 2) xtg_speak(sbn, 3, "Entering routine %s", sbn); |
| 166 | + |
| 167 | + zsam = (zmax - zmin) / (nzsam - 1); |
| 168 | + |
| 169 | + ib = 0; |
| 170 | + |
| 171 | + ibs1 = -1; |
| 172 | + ibs2 = -1; |
| 173 | + |
| 174 | + k1 = 1; |
| 175 | + k2 = nz; |
| 176 | + |
| 177 | + xtg_speak(sbn, 2, "Total number of XY poinst and Z points are %d %d", nxvec, nzsam); |
| 178 | + |
| 179 | + for (ic = 0; ic < nxvec; ic++) { |
| 180 | + xc = xvec[ic]; |
| 181 | + yc = yvec[ic]; |
| 182 | + if (debug > 1) xtg_speak(sbn, 2, "Column %d... X Y %f12.2 %f12.2", ic, xc, yc); |
| 183 | + |
| 184 | + _get_ij_range(&i1, &i2, &j1, &j2, xc, yc, mcol, mrow, xori, yori, xinc, yinc, |
| 185 | + yflip, rotation, maptopi, maptopj, mapbasi, mapbasj, debug); |
| 186 | + |
| 187 | + if (debug > 1) xtg_speak(sbn, 2, "I J range %d %d %d %d...", i1, i2, j1, j2); |
| 188 | + |
| 189 | + for (izc = 0; izc < nzsam; izc++) { |
| 190 | + |
| 191 | + zc = zmin + izc * zsam; |
| 192 | + |
| 193 | + /* check the onelayer version of the grid first (speed up) */ |
| 194 | + ier = grd3d_point_val_crange(xc, yc, zc, nx, ny, 1, p_coor_v, |
| 195 | + p_zcornone_v, p_actnumone_v, p_dummy_v, &value, |
| 196 | + i1, i2, j1, j2, 1, 1, &ibs1, -1, debug); |
| 197 | + |
| 198 | + |
| 199 | + if (ier == 0) { |
| 200 | + |
| 201 | + if (debug > 2 && ier == 0) xtg_speak(sbn, 3, "Trying K1 K2 %d %d", |
| 202 | + k1, k2); |
| 203 | + |
| 204 | + ios = grd3d_point_val_crange(xc, yc, zc, nx, ny, nz, p_coor_v, |
| 205 | + p_zcorn_v, p_actnum_v, p_val_v, |
| 206 | + &value, i1, i2, j1, j2, k1, k2, |
| 207 | + &ibs2, 0, debug); |
| 208 | + |
| 209 | + if (ios == 0) { |
| 210 | + values[ib++] = value; |
| 211 | + } |
| 212 | + else{ |
| 213 | + values[ib++] = UNDEF; |
| 214 | + } |
| 215 | + |
| 216 | + } |
| 217 | + else{ |
| 218 | + /* outside onelayer cell */ |
| 219 | + values[ib++] = UNDEF; |
| 220 | + continue; |
| 221 | + } |
| 222 | + } |
| 223 | + } |
| 224 | + |
| 225 | + if (debug > 2) xtg_speak(sbn, 3, "Done ..."); |
| 226 | + |
| 227 | + return EXIT_SUCCESS; |
| 228 | + |
| 229 | +} |
0 commit comments