1 /*************************************************************************** 2 triangulation.cpp - GDL library function 3 ------------------- 4 begin : Aug 30 2017 5 copyright : (C) 2017 by Gilles Duvert 6 7 ***************************************************************************/ 8 9 /*************************************************************************** 10 * * 11 * This program is free software; you can redistribute it and/or modify * 12 * it under the terms of the GNU General Public License as published by * 13 * the Free Software Foundation; either version 2 of the License, or * 14 * (at your option) any later version. * 15 * * 16 ***************************************************************************/ 17 18 #include "includefirst.hpp" 19 #include "datatypes.hpp" 20 #include "envt.hpp" 21 #include "dinterpreter.hpp" 22 23 using namespace std; 24 #include "delaunator/delaunator-header-only.hpp" 25 #include "tripack.c" 26 #include "stripack.c" 27 #include "ssrfpack.c" 28 //#include "akima760.c" 29 //#include "akima761.c" 30 31 namespace lib { 32 33 // using namespace std; 34 nextHalfedge(size_t e)35 inline size_t nextHalfedge(size_t e) { return (e % 3 == 2) ? e - 2 : e + 1; } prevHalfedge(size_t e)36 inline size_t prevHalfedge(size_t e) { return (e % 3 == 0) ? e + 2 : e - 1; } triangleOfEdge(size_t e)37 inline size_t triangleOfEdge(size_t e) { return floor(e / 3); } 38 GDL_Triangulate(EnvT * e)39 void GDL_Triangulate(EnvT* e) { 40 41 DDoubleGDL *xVal, *yVal, *fvalue; 42 DLong npts; 43 SizeT nParam = e->NParam(); 44 if (nParam < 2) e->Throw("Incorrect number of arguments."); //actually IDL permits to *not* have the 3rd argument. 45 bool wantsTriangles=(nParam > 2); 46 bool wantsEdge=(nParam == 4); 47 if (wantsTriangles) e->AssureGlobalPar(2); //since we return values in it? 48 if (wantsEdge) e->AssureGlobalPar(3); //since we return values in it? 49 50 static int sphereIx=e->KeywordIx( "SPHERE"); 51 static int degreeIx=e->KeywordIx( "DEGREES"); 52 static int fvalueIx=e->KeywordIx( "FVALUE"); 53 bool isSphere=(e->KeywordPresent(sphereIx)); 54 bool doDegree=(e->KeywordPresent(degreeIx)); 55 bool hasFvalue=(e->KeywordPresent(fvalueIx)); 56 if (!isSphere) { 57 doDegree=false; 58 hasFvalue=false; 59 } 60 bool wantsDupes=false; 61 static int dupesIx=e->KeywordIx( "REPEATS"); 62 if( e->KeywordPresent( dupesIx)) wantsDupes=true; 63 64 bool wantsConnectivity=false; 65 static int connIx=e->KeywordIx( "CONNECTIVITY"); 66 if( e->KeywordPresent( connIx)) wantsConnectivity=true; 67 68 //check since we return values in them 69 if (isSphere) e->AssureGlobalKW(sphereIx); 70 if (hasFvalue) e->AssureGlobalKW(fvalueIx); 71 if (wantsDupes) e->AssureGlobalKW(dupesIx); 72 if (wantsConnectivity) e->AssureGlobalKW(connIx); 73 74 //some things depend on X and Y being doubles of not 75 bool isDouble=false; 76 if ( e->GetParDefined(0)->Type()==GDL_DOUBLE && e->GetParDefined(1)->Type()==GDL_DOUBLE ) isDouble=true; 77 xVal = e->GetParAs< DDoubleGDL > (0); 78 if (xVal->Rank() == 0) e->Throw("Expression must be an array in this context: " + e->GetParString(0)); 79 npts = xVal->N_Elements(); 80 yVal = e->GetParAs< DDoubleGDL > (1); 81 if (yVal->Rank() == 0) e->Throw("Expression must be an array in this context: " + e->GetParString(1)); 82 if (yVal->N_Elements() != npts) e->Throw("X & Y arrays must have same number of points."); 83 if (hasFvalue) { 84 fvalue=e->GetKWAs<DDoubleGDL>(fvalueIx); 85 if (fvalue->Rank() == 0) e->Throw("Expression must be an array in this context: " + e->GetString( fvalueIx)+"."); 86 if (fvalue->N_Elements() != npts) e->Throw("X & Y arrays must have same number of points."); //yes yes. 87 } 88 89 DLong maxEl, minEl; 90 DDouble maxVal; 91 xVal->MinMax(&minEl,&maxEl,NULL,NULL,false); 92 //maximum ABSOLUTE value 93 maxVal=abs((*xVal)[maxEl]); 94 yVal->MinMax(&minEl,&maxEl,NULL,NULL,false); 95 //maximum ABSOLUTE value 96 maxVal=std::max(maxVal,abs((*yVal)[maxEl])); 97 98 DDouble dtol = isDouble ? 1e-12 : 1e-6; 99 //Tol is irrelevant in our implementation, as (s)tripack work with tol=machine precision. 100 101 DDouble* xx=&(*xVal)[0]; 102 DDouble* yy=&(*yVal)[0]; 103 104 //IN SPHERE MODE, xVal and yVal ARE RETURNED, ARE DOUBLE PRECISION and their ORDER is MODIFIED, 105 //the input points are sorted using the coordinate (x or y) covering max range (y prefered). It is not the case here. 106 //SPHERE does not support duplicate points yet. 107 108 if (isSphere) { 109 DDoubleGDL *sc_xVal=xVal; 110 DDoubleGDL *sc_yVal=yVal; 111 Guard<BaseGDL>xguard; 112 Guard<BaseGDL>yguard; 113 if (doDegree) { 114 static DDouble DToR=double(3.1415926535897932384626433832795)/180.0; 115 sc_xVal=new DDoubleGDL(npts,BaseGDL::NOZERO); 116 for (DLong i=0; i<npts; ++i) (*sc_xVal)[i]=(*xVal)[i]*DToR; 117 xguard.reset(sc_xVal); 118 sc_yVal=new DDoubleGDL(npts,BaseGDL::NOZERO); 119 for (DLong i=0; i<npts; ++i) (*sc_yVal)[i]=(*yVal)[i]*DToR; 120 yguard.reset(sc_yVal); 121 } 122 DDoubleGDL* x=new DDoubleGDL(npts,BaseGDL::NOZERO); 123 DDoubleGDL* y=new DDoubleGDL(npts,BaseGDL::NOZERO); 124 DDoubleGDL* z=new DDoubleGDL(npts,BaseGDL::NOZERO); 125 //Stripack is lat, lon and IDL lon,lat 126 DLong ret0=stripack::trans_(&npts, (DDouble*)sc_yVal->DataAddr(), (DDouble*)sc_xVal->DataAddr(),(DDouble*)x->DataAddr(),(DDouble*)y->DataAddr(),(DDouble*)z->DataAddr()); 127 DLong listsize=6*npts-12; 128 //we use GDL objects as they will be visible outside at GDL level. 129 DLongGDL* list=new DLongGDL(listsize,BaseGDL::NOZERO); 130 DLongGDL* lptr=new DLongGDL(listsize,BaseGDL::NOZERO); 131 DLongGDL* lend=new DLongGDL(npts,BaseGDL::NOZERO); 132 DLong* near__=(DLong*)malloc(npts*sizeof(DLong)); //initial name "near" would not work as this is reserved in Windows C. 133 DLong* next=(DLong*)malloc(npts*sizeof(DLong)); 134 DDouble* dist=(DDouble*)malloc(npts*sizeof(DDouble)); 135 DLong ier=0; 136 DLong lnew=0; 137 DLong ret1=stripack::sph_trmesh_(&npts,(DDouble*)x->DataAddr(),(DDouble*)y->DataAddr(),(DDouble*)z->DataAddr(), (DLong*)list->DataAddr(), (DLong*)lptr->DataAddr(), (DLong*)lend->DataAddr(), &lnew, near__, next, dist, &ier); 138 free(dist); 139 free(next); 140 free(near__); 141 if (ier !=0) { 142 GDLDelete(list); 143 GDLDelete(lptr); 144 GDLDelete(lend); 145 GDLDelete(x); 146 GDLDelete(y); 147 GDLDelete(z); 148 if (ier == -2) e->Throw("The first 3 nodes lie on the same great circle."); 149 else e->Throw("Duplicate nodes encountered."); 150 } 151 //Structure returned in SPHERE is not the same as IDL (although it must be possible to write 152 //something identical, but for what use? The structure will be passed to the TRIGRID function 153 //that uses SRFPACK, which in turn uses this structure, containing x,y,z,list,lptr,lend 154 // creating the output anonymous structure 155 DStructDesc* sphere_desc = new DStructDesc("$truct"); 156 SpDDouble aDbleArr(dimension((SizeT)npts)); 157 SpDLong aLongArr1(dimension((SizeT)listsize)); 158 SpDLong aLongArr2(dimension((SizeT)npts)); 159 sphere_desc->AddTag("X", &aDbleArr); 160 sphere_desc->AddTag("Y", &aDbleArr); 161 sphere_desc->AddTag("Z", &aDbleArr); 162 sphere_desc->AddTag("LIST", &aLongArr1); 163 sphere_desc->AddTag("LPTR", &aLongArr1); 164 sphere_desc->AddTag("LEND", &aLongArr2); 165 DStructGDL* sphere = new DStructGDL(sphere_desc, dimension()); 166 167 sphere->InitTag("X", *x); 168 sphere->InitTag("Y", *y); 169 sphere->InitTag("Z", *z); 170 sphere->InitTag("LIST", *list); 171 sphere->InitTag("LPTR", *lptr); 172 sphere->InitTag("LEND", *lend); 173 e->SetKW(sphereIx,sphere); 174 175 if (wantsTriangles) { 176 //convert to triangle list 177 DLong nrow = 6; //no arcs. 178 DLong* ltri = (DLong*) malloc((12 * npts) * sizeof (DLong)); 179 DLong ntriangles = 0; 180 DLong ret2=stripack::sph_trlist_(&npts, (DLong*)list->DataAddr(), (DLong*)lptr->DataAddr(), (DLong*)lend->DataAddr(), &nrow, &ntriangles, ltri, &ier); 181 if (ier != 0) 182 { 183 free(ltri); 184 e->Throw("Unexpected Error in STRIPACK, SPH_TRLIST routine. Please report."); 185 } 186 SizeT d[2]; 187 d[1] = ntriangles; 188 d[0] = 3; 189 DLongGDL* returned_triangles = new DLongGDL(dimension(d, 2), BaseGDL::NOZERO); 190 for (SizeT j = 0; j < ntriangles; ++j) 191 { 192 for (int i = 0; i < 3; ++i) 193 { 194 (*returned_triangles)[3 * j + i] = (ltri[6 * j + i] - 1); //nubering starts a 0. 195 } 196 } 197 free(ltri); 198 //pass back to GDL env: 199 e->SetPar(2, returned_triangles); 200 } 201 if (wantsEdge) { 202 DLong* nodes=(DLong*)malloc(npts*sizeof(DLong)); 203 DLong nb=0; 204 DLong na=0; 205 DLong nt=0; 206 DLong ret3=stripack::sph_bnodes_(&npts, (DLong*)list->DataAddr(), (DLong*)lptr->DataAddr(), (DLong*)lend->DataAddr(), nodes, &nb, &na, &nt); 207 if (nb>0){ 208 DLongGDL* returned_edges = new DLongGDL(nb, BaseGDL::NOZERO); 209 for (SizeT j = 0; j < nb; ++j) (*returned_edges)[j]=nodes[j]-1; 210 e->SetPar(3, returned_edges); 211 } else e->SetPar(3, new DLongGDL(-1)); 212 free(nodes); 213 } 214 if (hasFvalue) { 215 //No need to reorder whatever. 216 // //create a dummy array of same size 217 // DDoubleGDL* ret=new DDoubleGDL(npts,BaseGDL::NOZERO); 218 // //do not forget remove 1 to get C array indexes 219 // DLong index_lend, index_list; //can be negative. 220 // (*ret)[0]=(*fvalue)[(*list)[0]-1]; 221 // for (SizeT i=0; i<npts-1; ++i) { 222 // index_lend=(*lend)[i]; //as it starts at 1, is the next point, not the end of current point neighbours. 223 // index_list=(*list)[index_lend]; //can be negative 224 // if (index_list < 0) index_list=-index_list-1; else index_list--; //remove negative values and set as C index. 225 // (*ret)[i+1]=(*fvalue)[index_list]; //next point 226 // } 227 // e->SetKW(fvalueIx,ret); 228 } 229 //no connectivity in SPHERE mode in IDL, but yes we can! 230 if (wantsConnectivity) { 231 //remove 1 to get C array indexes. 232 for (SizeT i = 0; i < lnew-1; ++i) (*lptr)[i]--; 233 for (SizeT i = 0; i < npts; ++i) (*lend)[i]--; 234 DLong* array=(DLong*)malloc((2*(lnew-1))*sizeof(DLong)); // size > max possible connectivity 235 SizeT runningindex=npts+1; 236 SizeT startindex=0; 237 array[startindex++]=npts+1; 238 for (SizeT i = 0; i < npts; ++i) { 239 //is it an exterior point? Yes if the termination of the connectivity list is exterior. 240 DLong lpl=(*lend)[i]; 241 if ((*list)[lpl]<0) array[runningindex++]=i; //exterior - we write it 242 //write all points until nfin=lend[i] is found again using lptr connectivity pointers: 243 DLong lp=lpl; //see nbcnt_() 244 do { 245 lp=(*lptr)[lp]; 246 array[runningindex++]=((*list)[lp]>0)?(*list)[lp]-1:(-(*list)[lp])-1; 247 } while (lp!=lpl); 248 array[startindex++]=runningindex; 249 } 250 DLongGDL* connections = new DLongGDL(runningindex, BaseGDL::NOZERO); 251 for (SizeT i = 0; i < runningindex; ++i) (*connections)[i]=array[i]; 252 e->SetKW(connIx,connections); 253 free(array); 254 } 255 //no more cleanup, x,y,z,and list,lptr,lend are in the returned structure! 256 } else { 257 258 // // for PLANE triangulation, everything must be scaled in order to have triangulation independent of range, 259 // // as the triangulation code IS sensitive to the values of X and Y. 260 // if (maxVal > 0) { 261 // for (DLong i = 0; i < npts; ++i) xx[i] /= maxVal; 262 // for (DLong i = 0; i < npts; ++i) yy[i] /= maxVal; 263 // } 264 265 // for PLANE triangulation, it is possible that everything must be scaled in order to have triangulation independent of range. 266 // unfortunately this may have dire consequences. At the meoment I pass unscaled X and Y 267 std::vector<double> coords; 268 // for (DLong i = 0; i < npts; ++i) {coords.push_back(xx[i]/maxVal);coords.push_back(yy[i]/maxVal);} 269 for (DLong i = 0; i < npts; ++i) {coords.push_back(xx[i]);coords.push_back(yy[i]);} 270 delaunator::Delaunator tri(coords); 271 272 DLong ndupes=tri.dupes.size(); 273 if (wantsDupes) { 274 if (ndupes ==0) { 275 DLongGDL* nothing=new DLongGDL(dimension(2),BaseGDL::ZERO); 276 nothing->Dec(); 277 e->SetKW(dupesIx, nothing); 278 } 279 else { 280 DLongGDL* returned_dupes = new DLongGDL(dimension(2,ndupes), BaseGDL::NOZERO); 281 for (DLong i = 0; i < ndupes; ++i) { 282 (*returned_dupes)[2*i+0]=tri.dupes[i].first; 283 (*returned_dupes)[2*i+1]=tri.dupes[i].second; 284 } 285 e->SetKW(dupesIx, returned_dupes); 286 } 287 } 288 289 if (wantsTriangles) { 290 DLong ntriangles=tri.triangles.size()/3; 291 //convert to triangle list 292 SizeT d[2]; 293 d[1] = ntriangles; 294 d[0] = 3; 295 DLongGDL* returned_triangles = new DLongGDL(dimension(d, 2), BaseGDL::NOZERO); 296 for (DLong j = 0; j < ntriangles; ++j) 297 { 298 (*returned_triangles)[3 * j + 0] = tri.triangles[3*j+0]; 299 (*returned_triangles)[3 * j + 2] = tri.triangles[3*j+1]; 300 (*returned_triangles)[3 * j + 1] = tri.triangles[3*j+2]; //write it counterclockwise as IDL does. 301 } 302 e->SetPar(2, returned_triangles); 303 } 304 // get Hull values needed at two places. 305 std::vector<size_t> hull; 306 if (wantsEdge || wantsConnectivity) { 307 size_t index = tri.hull_start; 308 hull.push_back(index); 309 index = tri.hull_next[index]; 310 while (index != tri.hull_start) { 311 hull.push_back(index); 312 index = tri.hull_next[index]; 313 }; 314 } 315 if (wantsEdge) { 316 DLong nb=hull.size(); 317 DLongGDL* returned_edges = new DLongGDL(nb, BaseGDL::NOZERO); 318 for (DLong j = 0; j < nb; ++j) (*returned_edges)[j]=hull[nb-1-j]; //write it counterclockwise as IDL does. 319 e->SetPar(3, returned_edges); 320 } 321 if (wantsConnectivity) { 322 const std::size_t INVALID_INDEX = (std::numeric_limits<std::size_t>::max)(); 323 // as described in https://mapbox.github.io/delaunator/ it is not straightforward as there is no direct link points<->half-edges. 324 // let's build a map with a half-edge for each point, which will be a starting half-edge for connectivity. 325 //*** NOTE: CONNECTIVITY is oriented as in IDL 'naturally' *** 326 std::map<size_t,size_t>accel; //key is point, value is the out half-edge. 327 for (size_t e=0; e<tri.triangles.size();++e) { 328 size_t endpoint = tri.triangles[nextHalfedge(e)]; 329 if (accel.find(endpoint)==accel.end() || tri.halfedges[e]==INVALID_INDEX ) accel.insert(std::pair<size_t,size_t>(endpoint,e)); 330 } 331 // Now, points on hull will need a special treatment, as we must start the connectivity exploration at the previous hull position 332 // otherwise the connectivity search will end prematurely on a INVALID_INDEX 333 // We have computed the hull point list: for each hull point, find the halfedge that goes from the previous hull point to it. 334 //This is the halfedge to start with in accel: 335 for (int i=0; i< hull.size(); ++i) { 336 size_t point=hull[i]; //this is an existing point by construction, not a duplicate in dupes. 337 size_t nextpt=(i==0)?hull[hull.size()-1]:hull[i-1]; //idem. 338 std::map<size_t, size_t>::iterator iter = accel.find(nextpt); 339 size_t start = (*iter).second; //one of the halfedges starting with current point 340 size_t incoming=start; 341 //get adjacent halfedges until point is encountered. 342 size_t outgoing; 343 do { 344 outgoing = nextHalfedge(incoming); //the ougoing halfegde -> the point k in fact: assert(tri.triangles[outgoing]==effective_index[k]); 345 incoming = tri.halfedges[outgoing]; //the following incoming half-edge 346 } while (incoming!=INVALID_INDEX); //stop when halfedge going to 'point' is found. 347 iter = accel.find(point); 348 accel.erase(iter); 349 accel.insert(std::pair<size_t,size_t>(point,outgoing)); 350 } 351 352 //in connectivity we MUST have all the points, even the duplicated ones. Connectivity of duplicated points is wrong with IDL (intentional?). 353 // we could avoid this easily as I believe this is an IDL bug. We just have to reproduce the connectivity of the first encounter 354 // of the duplicated point. 355 DLong* array=(DLong*)malloc((npts*npts+npts+1)*sizeof(DLong)); // size > max possible connectivity 356 DLong runningindex = npts+1; // we report for all points 357 DLong startindex=0; 358 array[startindex++]= npts+1; 359 DLong* effective_index=(DLong*)malloc(npts*sizeof(DLong));//this is the list of npts vertexes for which we want the neighbour. It must be 360 // the list of l_npts (returned) indexes, with indexes of the first instance of duplicated points for the duplicated points. 361 for (DLong i=0; i< npts; ++i) effective_index[i]=i; 362 if (ndupes != 0) { 363 for (DLong idup=0; idup< ndupes; ++idup) { 364 DLong encounter=tri.dupes[idup].second; 365 effective_index[encounter]=tri.dupes[idup].first; 366 } 367 } 368 369 for (DLong k = 0; k < npts; ++k) { 370 //find the triangle in accel: 371 size_t point=effective_index[k]; 372 373 //if point is in the hull, make provision to add point itself at start and next point at end. 374 size_t nextpt; 375 bool isHull=false; 376 for (int i = 0; i < hull.size(); ++i) { 377 if (hull[i] == point) { 378 isHull=true; 379 nextpt = ((i + 1) >= hull.size()) ? hull[0] : hull[i + 1]; 380 array[runningindex++] = point; 381 break; 382 } 383 } 384 std::map<size_t, size_t>::iterator iter = accel.find(point); 385 size_t start = (*iter).second; //one of the halfedges starting with current point 386 size_t incoming=start; 387 //get adjacent halfedges and write their end points to connectivity until all encountered. 388 do { 389 /*THE FOLLOWING LINE IS KNOWN TO HAVE CRASHED ON UNREPRODUCTIBLE CIRCUMSTANCES ??? */ 390 array[runningindex++]=tri.triangles[incoming]; //the point outgoing halfedge ends. 391 size_t outgoing = nextHalfedge(incoming); //the ougoing halfegde -> the point k in fact: assert(tri.triangles[outgoing]==effective_index[k]); 392 incoming = tri.halfedges[outgoing]; //the following incoming half-edge 393 } while (incoming!=INVALID_INDEX && incoming != start); //stop when start halfedge is found. 394 if (isHull) array[runningindex++] = nextpt; 395 array[startindex++]=runningindex; 396 } 397 DLongGDL* connections = new DLongGDL(runningindex, BaseGDL::NOZERO); 398 for (DLong i = 0; i < runningindex; ++i) (*connections)[i]=array[i]; 399 e->SetKW(connIx,connections); 400 free(array); 401 } 402 } 403 } 404 405 trigrid_fun_spherical(EnvT * e)406 BaseGDL* trigrid_fun_spherical(EnvT* e) { 407 //NOT USED in this case: EXTRAPOLATE, MAX_VALUE, MIN_VALUE, QUINTIC, XOUT, YOUT 408 static DDouble DToR=double(3.1415926535897932384626433832795)/180.0; 409 SizeT nParam = e->NParam(); 410 if (nParam < 3) e->Throw("Incorrect number of arguments."); 411 //OK, trigrid does not care if more than 3 args in sphere mode. Limit at 6 is done at interpreter level. 412 413 // Get NX, NY values if present 414 DLong nx = 51; 415 DLong ny = 51; 416 417 static int degreeIx=e->KeywordIx( "DEGREES"); 418 bool doDegree=(e->KeywordPresent(degreeIx)); 419 420 // in difference with IDL we accept xout and yout for spherical data. 421 static int xoutIx=e->KeywordIx( "XOUT"); 422 bool doXout=(e->KeywordSet(xoutIx)); //and not "present" 423 static int youtIx=e->KeywordIx( "YOUT"); 424 bool doYout=(e->KeywordSet(youtIx)); 425 if ((doXout && !doYout) || (doYout && !doXout)) e->Throw("Incorrect number of arguments."); 426 static int xgridIx=e->KeywordIx( "XGRID"); 427 bool doXgrid=(e->KeywordPresent(xgridIx)); 428 static int ygridIx=e->KeywordIx( "YGRID"); 429 bool doYgrid=(e->KeywordPresent(ygridIx)); 430 431 static int nxIx = e->KeywordIx("NX"); 432 static int nyIx = e->KeywordIx("NY"); 433 bool canUseLimitsx=true; 434 bool canUseLimitsy=true; 435 if (e->KeywordSet(nxIx)) { canUseLimitsx=false; e->AssureLongScalarKW(nxIx, nx);} 436 if (e->KeywordSet(nyIx)) { canUseLimitsy=false; e->AssureLongScalarKW(nyIx, ny);} 437 if (nx < 1 || ny < 1) e->Throw("Array dimensions must be greater than 0."); 438 // we will further on define GS, whose value may overwrite the nx, ny. 439 440 BaseGDL* p0 = e->GetParDefined(0); //F 441 if (p0->Rank() == 0) e->Throw("Expression must be an array in this context: " + e->GetParString(0)); 442 DDoubleGDL* fval=static_cast<DDoubleGDL*>(p0->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 443 BaseGDL* p1 = e->GetParDefined(1); 444 if (p1->Rank() == 0) e->Throw("Expression must be an array in this context: " + e->GetParString(1)); 445 if (p1->N_Elements() != 2) e->Throw("Array must have 2 elements: " + e->GetParString(1)); 446 DDoubleGDL* GS = static_cast<DDoubleGDL*>(p1->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 447 448 BaseGDL* p2 = e->GetParDefined(2); 449 if (p2->Rank() == 0) e->Throw("Expression must be an array in this context: " + e->GetParString(2)); 450 if (p2->N_Elements() != 4) e->Throw("Array must have 4 elements: " + e->GetParString(2)); 451 DDoubleGDL* limits = static_cast<DDoubleGDL*> (p2->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 452 453 //Sphere is present since we are called. Test it. 454 DStructGDL* sphere=NULL; 455 static int sphereIx = e->KeywordIx( "SPHERE"); 456 BaseGDL* test=e->GetKW(sphereIx); 457 458 int xTag,yTag,zTag,listTag,lptrTag,lendTag; 459 460 if (test->Type() == GDL_STRUCT) 461 { 462 sphere = static_cast<DStructGDL*> (test); 463 464 xTag = sphere->Desc()->TagIndex("X"); 465 yTag = sphere->Desc()->TagIndex("Y"); 466 zTag = sphere->Desc()->TagIndex("Z"); 467 listTag = sphere->Desc()->TagIndex("LIST"); 468 lptrTag = sphere->Desc()->TagIndex("LPTR"); 469 lendTag = sphere->Desc()->TagIndex("LEND"); 470 if (xTag < 0 || yTag < 0 || zTag < 0 || listTag < 0 ||lptrTag < 0 || lendTag < 0) e->Throw("Invalid structure for the SPHERE keyword."); 471 } else { 472 e->Throw("SPHERE keyword must be a structure."); 473 } 474 475 DDoubleGDL* xVal = static_cast<DDoubleGDL*>(sphere->GetTag(xTag,0)); 476 DDoubleGDL* yVal = static_cast<DDoubleGDL*>(sphere->GetTag(yTag,0)); 477 DDoubleGDL* zVal = static_cast<DDoubleGDL*>(sphere->GetTag(zTag,0)); 478 if (xVal->N_Elements() != yVal->N_Elements())e->Throw("Invalid structure for the SPHERE keyword."); 479 if (xVal->N_Elements() != zVal->N_Elements())e->Throw("Invalid structure for the SPHERE keyword."); 480 DLongGDL* list = static_cast<DLongGDL*>(sphere->GetTag(listTag,0)); 481 DLongGDL* lptr = static_cast<DLongGDL*>(sphere->GetTag(lptrTag,0)); 482 DLongGDL* lend = static_cast<DLongGDL*>(sphere->GetTag(lendTag,0)); 483 484 //get npts 485 DLong npts=fval->N_Elements(); 486 if (xVal->N_Elements()<npts) npts=xVal->N_Elements(); 487 488 // Determine grid range 489 DDouble xref=0.0, yref=0.0; 490 DDouble xval, xinc; 491 DDouble yval, yinc; 492 DDouble xrange; 493 DDouble yrange; 494 //compute World positions of each pixels. 495 if ((*limits)[0]==(*limits)[2] || (*limits)[0]==(*limits)[2]) e->Throw("Inconsistent coordinate bounds."); 496 xval = (*limits)[0]; 497 if (doDegree) { 498 xval*=DToR; 499 xrange = (*limits)[2]*DToR - xval; 500 } else xrange = (*limits)[2] - xval; 501 yval = (*limits)[1]; 502 if (doDegree) { 503 yval*=DToR; 504 yrange = (*limits)[3]*DToR - yval; 505 } else yrange = (*limits)[3] - yval; 506 // Determine grid spacing 507 xinc = xrange / (nx-1); 508 yinc = yrange / (ny-1); 509 if (canUseLimitsx && (*GS)[0] > 0.0) 510 { 511 xinc = (*GS)[0]; 512 if (doDegree) xinc*=DToR; 513 nx = (xrange / xinc) +1; 514 } 515 if (canUseLimitsy && (*GS)[1] > 0.0) 516 { 517 yinc = (*GS)[1]; 518 if (doDegree) yinc*=DToR; 519 ny = (yrange / yinc) +1; 520 } 521 DDouble *plon, *plat; 522 bool cleanplon=false; 523 bool cleanplat=false; 524 525 if (doXout) { 526 DDoubleGDL* xOut=e->GetKWAs<DDoubleGDL>(xoutIx); 527 nx=xOut->N_Elements(); 528 if (nx < 1) e->Throw("XOUT must have >1 elements."); 529 plon=(DDouble*)xOut->DataAddr(); 530 if (doDegree) for (int i=0; i<nx; ++i) plon[i]*=DToR; 531 } else { 532 plon = (DDouble*)malloc(nx*sizeof(DDouble)); 533 for (int i=0; i<nx; ++i) plon[i]=(i-xref)*xinc+xval; 534 cleanplon=true; 535 } 536 537 if (doYout) { 538 DDoubleGDL* yOut=e->GetKWAs<DDoubleGDL>(youtIx); 539 ny=yOut->N_Elements(); 540 if (ny < 1) e->Throw("YOUT must have >1 elements."); 541 plat=(DDouble*)yOut->DataAddr(); 542 if (doDegree) for (int j=0; j<ny; ++j) plat[j]*=DToR; 543 } else { 544 plat = (DDouble*)malloc(ny*sizeof(DDouble)); 545 for (int j=0; j<ny; ++j) plat[j]=(j-yref)*yinc+yval; 546 cleanplat=true; 547 } 548 549 // Setup return array 550 DLong dims[2]; 551 dims[0] = ny; 552 dims[1] = nx; 553 dimension dim((DLong *) dims, 2); 554 DDoubleGDL* res = new DDoubleGDL(dim, BaseGDL::ZERO); 555 DLong ier=0; 556 557 //must be done here, unif_() below messes with plat and plon... 558 if (doXgrid) { 559 DDoubleGDL *xgrid=new DDoubleGDL(nx, BaseGDL::NOZERO); 560 for (int i=0; i<nx; ++i) (*xgrid)[i]=plon[i]; 561 if (doDegree) for (int i=0; i<nx; ++i) (*xgrid)[i]/=DToR; 562 e->SetKW(xgridIx, xgrid); 563 } 564 if (doYgrid) { 565 DDoubleGDL *ygrid=new DDoubleGDL(ny, BaseGDL::NOZERO); 566 for (int i=0; i<ny; ++i) (*ygrid)[i]=plat[i]; 567 if (doDegree) for (int i=0; i<ny; ++i) (*ygrid)[i]/=DToR; 568 e->SetKW(ygridIx, ygrid); 569 } 570 571 DLong iflgs=0; 572 DLong iflgg=0; 573 DDouble grad=0.0; 574 DDouble sigma=1.; 575 DLong ret2=ssrfpack::unif_ (&npts,(DDouble*)xVal->DataAddr(), (DDouble*)yVal->DataAddr(), 576 (DDouble*)zVal->DataAddr(), (DDouble*)fval->DataAddr(),(DLong*)list->DataAddr(), (DLong*)lptr->DataAddr(), 577 (DLong*)lend->DataAddr(), &iflgs, &sigma, &ny, &ny, &nx, plat, plon, &iflgg, 578 &grad, (DDouble*)res->DataAddr() , &ier); 579 if (cleanplon) free(plon); 580 if (cleanplat) free(plat); 581 if (ret2 != 0) { 582 if (ret2==-1) e->Throw("Error in SSRFPACK unif_() function: \"N, NI, NJ, or IFLGG is outside its valid range\"."); 583 if (ret2==-2) e->Throw("Error in SSRFPACK unif_() function: Nodes are colinear"); 584 if (ret2==-3) e->Throw("Error in SSRFPACK unif_() function: extrapolation failed due to the uniform grid extending too far beyond the triangulation boundary."); 585 } 586 587 return res->Transpose(NULL); 588 } 589 590 template< typename T1, typename T2> gdlGrid2DData(DLong nx,DDouble * x,DLong xref,DDouble xval,DDouble xinc,DLong ny,DDouble * y,DLong yref,DDouble yval,DDouble yinc,DLong ntri,DLongGDL * tri,DDoubleGDL * xVal,DDoubleGDL * yVal,T1 * zVal,bool domaxvalue,bool dominvalue,T2 maxVal,T2 minVal,T1 * res,bool input)591 void gdlGrid2DData(DLong nx, DDouble* x, DLong xref, DDouble xval, DDouble xinc, DLong ny, DDouble* y, DLong yref, DDouble yval, DDouble yinc, DLong ntri, DLongGDL* tri, 592 DDoubleGDL* xVal, DDoubleGDL* yVal, T1* zVal, bool domaxvalue, 593 bool dominvalue, T2 maxVal, T2 minVal, T1* res, bool input) { 594 595 if (dominvalue || domaxvalue) { 596 for (SizeT triIndex = 0; triIndex < ntri; ++triIndex) { 597 // first, find whose region of grid is concerned by this triangle: boundingBox. 598 DDouble xmin; 599 DDouble xmax; 600 DDouble ymin; 601 DDouble ymax; 602 DDouble x1, y1, x2, y2, x3, y3; 603 x1 = (*xVal)[(*tri)[3 * triIndex]]; 604 y1 = (*yVal)[(*tri)[3 * triIndex]]; 605 x2 = (*xVal)[(*tri)[3 * triIndex + 1]]; 606 y2 = (*yVal)[(*tri)[3 * triIndex + 1]]; 607 x3 = (*xVal)[(*tri)[3 * triIndex + 2]]; 608 y3 = (*yVal)[(*tri)[3 * triIndex + 2]]; 609 xmin = x1 < x2 ? (x1 < x3 ? x1 : x3) : (x2 < x3 ? x2 : x3); 610 ymin = y1 < y2 ? (y1 < y3 ? y1 : y3) : (y2 < y3 ? y2 : y3); 611 xmax = x1 > x2 ? (x1 > x3 ? x1 : x3) : (x2 > x3 ? x2 : x3); 612 ymax = y1 > y2 ? (y1 > y3 ? y1 : y3) : (y2 > y3 ? y2 : y3); 613 614 // x=(i-xref)*xinc+xval; -> i=(x-xval)/xinc+xref 615 int pxmin = ((xmin - xval) / xinc) + xref; 616 pxmin = (pxmin < 0) ? 0 : pxmin; 617 int pxmax = ((xmax - xval) / xinc) + xref + 1; 618 pxmax = (pxmax > nx) ? nx : pxmax; 619 pxmax = (pxmax < pxmin) ? pxmin : pxmax; 620 int pymin = ((ymin - yval) / yinc) + yref; 621 pymin = (pymin < 0) ? 0 : pymin; 622 int pymax = ((ymax - yval) / yinc) + yref + 1; 623 pymax = (pymax > ny) ? ny : pymax; 624 pymax = (pymax < pymin) ? pymin : pymax; 625 626 DDouble zc = (*zVal)[(*tri)[3 * triIndex + 2]]; 627 DDouble zac = (*zVal)[(*tri)[3 * triIndex]] - zc; 628 DDouble zbc = (*zVal)[(*tri)[3 * triIndex + 1]] - zc; 629 DDouble y23 = (y2 - y3); 630 DDouble x32 = (x3 - x2); 631 DDouble y31 = (y3 - y1); 632 DDouble x13 = (x1 - x3); 633 DDouble det = (y23 * x13 - x32 * y31); 634 DDouble minD = (det<0)?det:0 ; // we MUST be agnostic as to the orientation of the triangle, even if TRIANGULATE is clockwise. 635 DDouble maxD = (det>0)?det:0; 636 for (SizeT iy = pymin; iy < pymax; ++iy) { 637 DDouble dy = y[iy] - y3; 638 DLong start = iy*nx; 639 for (SizeT ix = pxmin; ix < pxmax; ++ix) { 640 DDouble dx = x[ix] - x3; 641 DDouble a = y23 * dx + x32 * dy; 642 if (a < minD || a > maxD) continue; 643 DDouble b = y31 * dx + x13 * dy; 644 if (b < minD || b > maxD) continue; 645 DDouble c = det - a - b; 646 if (c < minD || c > maxD) continue; 647 DDouble valbis = (a * zac + b * zbc) / det + zc; 648 if ((dominvalue && valbis < minVal) || (domaxvalue && valbis > maxVal)) { 649 } else (*res)[start + ix] = valbis; 650 } //ix loop 651 } // iy 652 } //tri 653 } else { 654 for (SizeT triIndex = 0; triIndex < ntri; ++triIndex) { 655 // first, find whose region of grid is concerned by this triangle: boundingBox. 656 DDouble xmin; 657 DDouble xmax; 658 DDouble ymin; 659 DDouble ymax; 660 DDouble x1, y1, x2, y2, x3, y3; 661 x1 = (*xVal)[(*tri)[3 * triIndex]]; 662 y1 = (*yVal)[(*tri)[3 * triIndex]]; 663 x2 = (*xVal)[(*tri)[3 * triIndex + 1]]; 664 y2 = (*yVal)[(*tri)[3 * triIndex + 1]]; 665 x3 = (*xVal)[(*tri)[3 * triIndex + 2]]; 666 y3 = (*yVal)[(*tri)[3 * triIndex + 2]]; 667 xmin = x1 < x2 ? (x1 < x3 ? x1 : x3) : (x2 < x3 ? x2 : x3); 668 ymin = y1 < y2 ? (y1 < y3 ? y1 : y3) : (y2 < y3 ? y2 : y3); 669 xmax = x1 > x2 ? (x1 > x3 ? x1 : x3) : (x2 > x3 ? x2 : x3); 670 ymax = y1 > y2 ? (y1 > y3 ? y1 : y3) : (y2 > y3 ? y2 : y3); 671 672 // x=(i-xref)*xinc+xval; -> i=(x-xval)/xinc+xref 673 int pxmin = ((xmin - xval) / xinc) + xref; 674 pxmin = (pxmin < 0) ? 0 : pxmin; 675 int pxmax = ((xmax - xval) / xinc) + xref + 1 ; 676 pxmax = (pxmax > nx ) ? nx : pxmax; 677 pxmax = (pxmax < pxmin) ? pxmin : pxmax; 678 int pymin = ((ymin - yval) / yinc) + yref; 679 pymin = (pymin < 0) ? 0 : pymin; 680 int pymax = ((ymax - yval) / yinc) + yref + 1; 681 pymax = (pymax > ny) ? ny : pymax; 682 pymax = (pymax < pymin) ? pymin : pymax; 683 684 DDouble zc = (*zVal)[(*tri)[3 * triIndex + 2]]; 685 DDouble zac = (*zVal)[(*tri)[3 * triIndex]] - zc; 686 DDouble zbc = (*zVal)[(*tri)[3 * triIndex + 1]] - zc; 687 DDouble y23 = (y2 - y3); 688 DDouble x32 = (x3 - x2); 689 DDouble y31 = (y3 - y1); 690 DDouble x13 = (x1 - x3); 691 DDouble det = (y23 * x13 - x32 * y31); 692 DDouble minD = (det<0)?det:0 ; // we MUST be agnostic as to the orientation of the triangle, even if TRIANGULATE is clockwise. 693 DDouble maxD = (det>0)?det:0; 694 for (SizeT iy = pymin; iy < pymax; ++iy) { 695 DDouble dy = y[iy]-y3; 696 DLong start = iy*nx; 697 for (SizeT ix = pxmin; ix < pxmax; ++ix) { 698 DDouble dx = x[ix]-x3; 699 DDouble a = y23 * dx + x32 * dy; 700 if (a < minD || a > maxD) continue; 701 DDouble b = y31 * dx + x13 * dy; 702 if (b < minD || b > maxD) continue; 703 DDouble c = det - a - b; 704 if (c < minD || c > maxD) continue; 705 (*res)[start + ix] = (a * zac + b * zbc)/det + zc; 706 } //ix loop 707 } // iy 708 } //tri 709 } 710 } 711 gdlGrid2DDataCpx(DLong nx,DDouble * x,DLong xref,DDouble xval,DDouble xinc,DLong ny,DDouble * y,DLong yref,DDouble yval,DDouble yinc,DLong ntri,DLongGDL * tri,DDoubleGDL * xVal,DDoubleGDL * yVal,DComplexDblGDL * zVal,bool domaxvalue,bool dominvalue,DComplexDbl maxVal,DComplexDbl minVal,DComplexDblGDL * res,bool input)712 void gdlGrid2DDataCpx(DLong nx, DDouble* x, DLong xref, DDouble xval, DDouble xinc, DLong ny, DDouble* y, DLong yref, DDouble yval, DDouble yinc, DLong ntri, DLongGDL* tri, 713 DDoubleGDL* xVal, DDoubleGDL* yVal, DComplexDblGDL* zVal, bool domaxvalue, 714 bool dominvalue, DComplexDbl maxVal, DComplexDbl minVal, DComplexDblGDL* res, bool input) { 715 716 if (dominvalue || domaxvalue) { 717 for (SizeT triIndex = 0; triIndex < ntri; ++triIndex) { 718 // first, find whose region of grid is concerned by this triangle: boundingBox. 719 DDouble xmin; 720 DDouble xmax; 721 DDouble ymin; 722 DDouble ymax; 723 DDouble x1, y1, x2, y2, x3, y3; 724 x1 = (*xVal)[(*tri)[3 * triIndex]]; 725 y1 = (*yVal)[(*tri)[3 * triIndex]]; 726 x2 = (*xVal)[(*tri)[3 * triIndex + 1]]; 727 y2 = (*yVal)[(*tri)[3 * triIndex + 1]]; 728 x3 = (*xVal)[(*tri)[3 * triIndex + 2]]; 729 y3 = (*yVal)[(*tri)[3 * triIndex + 2]]; 730 xmin = x1 < x2 ? (x1 < x3 ? x1 : x3) : (x2 < x3 ? x2 : x3); 731 ymin = y1 < y2 ? (y1 < y3 ? y1 : y3) : (y2 < y3 ? y2 : y3); 732 xmax = x1 > x2 ? (x1 > x3 ? x1 : x3) : (x2 > x3 ? x2 : x3); 733 ymax = y1 > y2 ? (y1 > y3 ? y1 : y3) : (y2 > y3 ? y2 : y3); 734 735 // x=(i-xref)*xinc+xval; -> i=(x-xval)/xinc+xref 736 int pxmin = ((xmin - xval) / xinc) + xref; 737 pxmin = (pxmin < 0) ? 0 : pxmin; 738 int pxmax = ((xmax - xval) / xinc) + xref + 1; 739 pxmax = (pxmax > nx) ? nx : pxmax; 740 pxmax = (pxmax < pxmin) ? pxmin : pxmax; 741 int pymin = ((ymin - yval) / yinc) + yref; 742 pymin = (pymin < 0) ? 0 : pymin; 743 int pymax = ((ymax - yval) / yinc) + yref + 1; 744 pymax = (pymax > ny) ? ny : pymax; 745 pymax = (pymax < pymin) ? pymin : pymax; 746 747 DComplexDbl zc = (*zVal)[(*tri)[3 * triIndex + 2]]; 748 DComplexDbl zac = std::complex<double>((*zVal)[(*tri)[3 * triIndex]].real() - zc.real(), (*zVal)[(*tri)[3 * triIndex]].imag() - zc.imag()); 749 DComplexDbl zbc = std::complex<double>((*zVal)[(*tri)[3 * triIndex + 1]].real() - zc.real(), (*zVal)[(*tri)[3 * triIndex + 1]].imag() - zc.imag()); 750 DDouble y23 = (y2 - y3); 751 DDouble x32 = (x3 - x2); 752 DDouble y31 = (y3 - y1); 753 DDouble x13 = (x1 - x3); 754 DDouble det = (y23 * x13 - x32 * y31); 755 DDouble minD = (det<0)?det:0 ; // we MUST be agnostic as to the orientation of the triangle, even if TRIANGULATE is clockwise. 756 DDouble maxD = (det>0)?det:0; 757 for (SizeT iy = pymin; iy < pymax; ++iy) { 758 DDouble dy = y[iy] - y3; 759 DLong start = iy*nx; 760 for (SizeT ix = pxmin; ix < pxmax; ++ix) { 761 DDouble dx = x[ix] - x3; 762 DDouble a = y23 * dx + x32 * dy; 763 if (a < minD || a > maxD) continue; 764 DDouble b = y31 * dx + x13 * dy; 765 if (b < minD || b > maxD) continue; 766 DDouble c = det - a - b; 767 if (c < minD || c > maxD) continue; 768 DDouble valbisR = (a * zac.real() + b * zbc.real())/det + zc.real(); 769 DDouble valbisI = (a * zac.imag() + b * zbc.imag())/det + zc.imag(); 770 if ((dominvalue && valbisR < minVal.real()) || (domaxvalue && valbisR > maxVal.real())) { valbisR = (*res)[start + ix].real(); } 771 if ((dominvalue && valbisI < minVal.imag()) || (domaxvalue && valbisI > maxVal.imag())) { valbisI = (*res)[start + ix].imag(); } 772 (*res)[start + ix] = std::complex<double>(valbisR,valbisI); 773 } //ix loop 774 } // iy 775 } //tri 776 } else { 777 for (SizeT triIndex = 0; triIndex < ntri; ++triIndex) { 778 // first, find whose region of grid is concerned by this triangle: boundingBox. 779 DDouble xmin; 780 DDouble xmax; 781 DDouble ymin; 782 DDouble ymax; 783 DDouble x1, y1, x2, y2, x3, y3; 784 x1 = (*xVal)[(*tri)[3 * triIndex]]; 785 y1 = (*yVal)[(*tri)[3 * triIndex]]; 786 x2 = (*xVal)[(*tri)[3 * triIndex + 1]]; 787 y2 = (*yVal)[(*tri)[3 * triIndex + 1]]; 788 x3 = (*xVal)[(*tri)[3 * triIndex + 2]]; 789 y3 = (*yVal)[(*tri)[3 * triIndex + 2]]; 790 xmin = x1 < x2 ? (x1 < x3 ? x1 : x3) : (x2 < x3 ? x2 : x3); 791 ymin = y1 < y2 ? (y1 < y3 ? y1 : y3) : (y2 < y3 ? y2 : y3); 792 xmax = x1 > x2 ? (x1 > x3 ? x1 : x3) : (x2 > x3 ? x2 : x3); 793 ymax = y1 > y2 ? (y1 > y3 ? y1 : y3) : (y2 > y3 ? y2 : y3); 794 795 // x=(i-xref)*xinc+xval; -> i=(x-xval)/xinc+xref 796 int pxmin = ((xmin - xval) / xinc) + xref; 797 pxmin = (pxmin < 0) ? 0 : pxmin; 798 int pxmax = ((xmax - xval) / xinc) + xref + 1; 799 pxmax = (pxmax > nx) ? nx : pxmax; 800 pxmax = (pxmax < pxmin) ? pxmin : pxmax; 801 int pymin = ((ymin - yval) / yinc) + yref; 802 pymin = (pymin < 0) ? 0 : pymin; 803 int pymax = ((ymax - yval) / yinc) + yref + 1; 804 pymax = (pymax > ny) ? ny : pymax; 805 pymax = (pymax < pymin) ? pymin : pymax; 806 807 DComplexDbl zc = (*zVal)[(*tri)[3 * triIndex + 2]]; 808 DComplexDbl zac = std::complex<double>((*zVal)[(*tri)[3 * triIndex]].real() - zc.real(), (*zVal)[(*tri)[3 * triIndex]].imag() - zc.imag()); 809 DComplexDbl zbc = std::complex<double>((*zVal)[(*tri)[3 * triIndex + 1]].real() - zc.real(), (*zVal)[(*tri)[3 * triIndex + 1]].imag() - zc.imag()); 810 DDouble y23 = (y2 - y3); 811 DDouble x32 = (x3 - x2); 812 DDouble y31 = (y3 - y1); 813 DDouble x13 = (x1 - x3); 814 DDouble det = (y23 * x13 - x32 * y31); 815 DDouble minD = (det<0)?det:0 ; // we MUST be agnostic as to the orientation of the triangle, even if TRIANGULATE is clockwise. 816 DDouble maxD = (det>0)?det:0; 817 for (SizeT iy = pymin; iy < pymax; ++iy) { 818 DDouble dy = y[iy] - y3; 819 DLong start = iy*nx; 820 for (SizeT ix = pxmin; ix < pxmax; ++ix) { 821 DDouble dx = x[ix] - x3; 822 DDouble a = y23 * dx + x32 * dy; 823 if (a < minD || a > maxD) continue; 824 DDouble b = y31 * dx + x13 * dy; 825 if (b < minD || b > maxD) continue; 826 DDouble c = det - a - b; 827 if (c < minD || c > maxD) continue; 828 DDouble valbisR = (a * zac.real() + b * zbc.real())/det + zc.real(); 829 DDouble valbisI = (a * zac.imag() + b * zbc.imag())/det + zc.imag(); 830 (*res)[start + ix] = std::complex<double>(valbisR,valbisI); 831 } //ix loop 832 } // iy 833 } //tri 834 } 835 } 836 837 template< typename T1, typename T2> gdlGrid2DData(DLong nx,DDouble * x,DLong ny,DDouble * y,DLong ntri,DLongGDL * tri,DDoubleGDL * xVal,DDoubleGDL * yVal,T1 * zVal,bool domaxvalue,bool dominvalue,T2 maxVal,T2 minVal,T1 * res,bool input)838 void gdlGrid2DData(DLong nx, DDouble* x, DLong ny, DDouble* y, DLong ntri, DLongGDL* tri, 839 DDoubleGDL* xVal, DDoubleGDL* yVal, T1* zVal, bool domaxvalue, 840 bool dominvalue, T2 maxVal, T2 minVal, T1* res, bool input){ 841 842 if (dominvalue || domaxvalue) { 843 for (SizeT triIndex = 0; triIndex < ntri; ++triIndex) { 844 // first, find whose region of grid is concerned by this triangle: boundingBox. 845 DDouble xmin; 846 DDouble xmax; 847 DDouble ymin; 848 DDouble ymax; 849 DDouble x1, y1, x2, y2, x3, y3; 850 x1 = (*xVal)[(*tri)[3 * triIndex]]; 851 y1 = (*yVal)[(*tri)[3 * triIndex]]; 852 x2 = (*xVal)[(*tri)[3 * triIndex + 1]]; 853 y2 = (*yVal)[(*tri)[3 * triIndex + 1]]; 854 x3 = (*xVal)[(*tri)[3 * triIndex + 2]]; 855 y3 = (*yVal)[(*tri)[3 * triIndex + 2]]; 856 xmin = x1 < x2 ? (x1 < x3 ? x1 : x3) : (x2 < x3 ? x2 : x3); 857 ymin = y1 < y2 ? (y1 < y3 ? y1 : y3) : (y2 < y3 ? y2 : y3); 858 xmax = x1 > x2 ? (x1 > x3 ? x1 : x3) : (x2 > x3 ? x2 : x3); 859 ymax = y1 > y2 ? (y1 > y3 ? y1 : y3) : (y2 > y3 ? y2 : y3); 860 861 // find first and last index of x and y that will intersect this BBox 862 DLong k = 0; 863 while (k < (nx - 1) && x[k + 1] < xmin) { 864 k++; 865 } 866 int pxmin = k; 867 while (k < (nx - 1) && x[k + 1] < xmax) { 868 k++; 869 } 870 int pxmax = k; 871 k = 0; 872 while (k < (ny - 1) && y[k + 1] < ymin) { 873 k++; 874 } 875 int pymin = k; 876 while (k < (ny - 1) && y[k + 1] < ymax) { 877 k++; 878 } 879 int pymax = k; 880 881 DDouble zc = (*zVal)[(*tri)[3 * triIndex + 2]]; 882 DDouble zac = (*zVal)[(*tri)[3 * triIndex]] - zc; 883 DDouble zbc = (*zVal)[(*tri)[3 * triIndex + 1]] - zc; 884 DDouble y23 = (y2 - y3); 885 DDouble x32 = (x3 - x2); 886 DDouble y31 = (y3 - y1); 887 DDouble x13 = (x1 - x3); 888 DDouble det = (y23 * x13 - x32 * y31); 889 DDouble minD = (det<0)?det:0 ; // we MUST be agnostic as to the orientation of the triangle, even if TRIANGULATE is clockwise. 890 DDouble maxD = (det>0)?det:0; 891 for (SizeT iy = pymin; iy <= pymax; ++iy) { 892 DDouble dy = y[iy] - y3; 893 DLong start = iy*nx; 894 for (SizeT ix = pxmin; ix <= pxmax; ++ix) { 895 DDouble dx = x[ix] - x3; 896 DDouble a = y23 * dx + x32 * dy; 897 if (a < minD || a > maxD) continue; 898 DDouble b = y31 * dx + x13 * dy; 899 if (b < minD || b > maxD) continue; 900 DDouble c = det - a - b; 901 if (c < minD || c > maxD) continue; 902 DDouble valbis = (a * zac + b * zbc) / det + zc; 903 if ((dominvalue && valbis < minVal) || (domaxvalue && valbis > maxVal)) { 904 } else (*res)[start + ix] = valbis; 905 } //ix loop 906 } // iy 907 } //tri 908 } else { 909 for (SizeT triIndex = 0; triIndex < ntri; ++triIndex) { 910 // first, find whose region of grid is concerned by this triangle: boundingBox. 911 DDouble xmin; 912 DDouble xmax; 913 DDouble ymin; 914 DDouble ymax; 915 DDouble x1, y1, x2, y2, x3, y3; 916 x1 = (*xVal)[(*tri)[3 * triIndex]]; 917 y1 = (*yVal)[(*tri)[3 * triIndex]]; 918 x2 = (*xVal)[(*tri)[3 * triIndex + 1]]; 919 y2 = (*yVal)[(*tri)[3 * triIndex + 1]]; 920 x3 = (*xVal)[(*tri)[3 * triIndex + 2]]; 921 y3 = (*yVal)[(*tri)[3 * triIndex + 2]]; 922 xmin = x1 < x2 ? (x1 < x3 ? x1 : x3) : (x2 < x3 ? x2 : x3); 923 ymin = y1 < y2 ? (y1 < y3 ? y1 : y3) : (y2 < y3 ? y2 : y3); 924 xmax = x1 > x2 ? (x1 > x3 ? x1 : x3) : (x2 > x3 ? x2 : x3); 925 ymax = y1 > y2 ? (y1 > y3 ? y1 : y3) : (y2 > y3 ? y2 : y3); 926 927 // find first and last index of x and y that will intersect this BBox 928 DLong k = 0; 929 while (k < (nx - 1) && x[k + 1] < xmin) { 930 k++; 931 } 932 int pxmin = k; 933 while (k < (nx - 1) && x[k + 1] < xmax) { 934 k++; 935 } 936 int pxmax = k; 937 k = 0; 938 while (k < (ny - 1) && y[k + 1] < ymin) { 939 k++; 940 } 941 int pymin = k; 942 while (k < (ny - 1) && y[k + 1] < ymax) { 943 k++; 944 } 945 int pymax = k; 946 947 DDouble zc = (*zVal)[(*tri)[3 * triIndex + 2]]; 948 DDouble zac = (*zVal)[(*tri)[3 * triIndex]] - zc; 949 DDouble zbc = (*zVal)[(*tri)[3 * triIndex + 1]] - zc; 950 DDouble y23 = (y2 - y3); 951 DDouble x32 = (x3 - x2); 952 DDouble y31 = (y3 - y1); 953 DDouble x13 = (x1 - x3); 954 DDouble det = (y23 * x13 - x32 * y31); 955 DDouble minD = (det<0)?det:0 ; // we MUST be agnostic as to the orientation of the triangle, even if TRIANGULATE is clockwise. 956 DDouble maxD = (det>0)?det:0; 957 for (SizeT iy = pymin; iy <= pymax; ++iy) { 958 DDouble dy = y[iy]-y3; 959 DLong start = iy*nx; 960 for (SizeT ix = pxmin; ix <= pxmax; ++ix) { 961 DDouble dx = x[ix]-x3; 962 DDouble a = y23 * dx + x32 * dy; 963 if (a < minD || a > maxD) continue; 964 DDouble b = y31 * dx + x13 * dy; 965 if (b < minD || b > maxD) continue; 966 DDouble c = det - a - b; 967 if (c < minD || c > maxD) continue; 968 (*res)[start + ix] = (a * zac + b * zbc)/det + zc; 969 } //ix loop 970 } // iy 971 } //tri 972 } 973 } 974 gdlGrid2DDataCpx(DLong nx,DDouble * x,DLong ny,DDouble * y,DLong ntri,DLongGDL * tri,DDoubleGDL * xVal,DDoubleGDL * yVal,DComplexDblGDL * zVal,bool domaxvalue,bool dominvalue,DComplexDbl maxVal,DComplexDbl minVal,DComplexDblGDL * res,bool input)975 void gdlGrid2DDataCpx(DLong nx, DDouble* x, DLong ny, DDouble* y, DLong ntri, DLongGDL* tri, 976 DDoubleGDL* xVal, DDoubleGDL* yVal, DComplexDblGDL* zVal, bool domaxvalue, 977 bool dominvalue, DComplexDbl maxVal, DComplexDbl minVal, DComplexDblGDL* res, bool input) { 978 979 if (dominvalue || domaxvalue) { 980 for (SizeT triIndex = 0; triIndex < ntri; ++triIndex) { 981 // first, find whose region of grid is concerned by this triangle: boundingBox. 982 DDouble xmin; 983 DDouble xmax; 984 DDouble ymin; 985 DDouble ymax; 986 DDouble x1, y1, x2, y2, x3, y3; 987 x1 = (*xVal)[(*tri)[3 * triIndex]]; 988 y1 = (*yVal)[(*tri)[3 * triIndex]]; 989 x2 = (*xVal)[(*tri)[3 * triIndex + 1]]; 990 y2 = (*yVal)[(*tri)[3 * triIndex + 1]]; 991 x3 = (*xVal)[(*tri)[3 * triIndex + 2]]; 992 y3 = (*yVal)[(*tri)[3 * triIndex + 2]]; 993 xmin = x1 < x2 ? (x1 < x3 ? x1 : x3) : (x2 < x3 ? x2 : x3); 994 ymin = y1 < y2 ? (y1 < y3 ? y1 : y3) : (y2 < y3 ? y2 : y3); 995 xmax = x1 > x2 ? (x1 > x3 ? x1 : x3) : (x2 > x3 ? x2 : x3); 996 ymax = y1 > y2 ? (y1 > y3 ? y1 : y3) : (y2 > y3 ? y2 : y3); 997 998 // find first and last index of x and y that will intersect this BBox 999 DLong k = 0; 1000 while (k < (nx - 1) && x[k + 1] < xmin) { 1001 k++; 1002 } 1003 int pxmin = k; 1004 while (k < (nx - 1) && x[k + 1] < xmax) { 1005 k++; 1006 } 1007 int pxmax = k; 1008 k = 0; 1009 while (k < (ny - 1) && y[k + 1] < ymin) { 1010 k++; 1011 } 1012 int pymin = k; 1013 while (k < (ny - 1) && y[k + 1] < ymax) { 1014 k++; 1015 } 1016 int pymax = k; 1017 1018 DComplexDbl zc = (*zVal)[(*tri)[3 * triIndex + 2]]; 1019 DComplexDbl zac = std::complex<double>((*zVal)[(*tri)[3 * triIndex]].real() - zc.real(), (*zVal)[(*tri)[3 * triIndex]].imag() - zc.imag()); 1020 DComplexDbl zbc = std::complex<double>((*zVal)[(*tri)[3 * triIndex + 1]].real() - zc.real(), (*zVal)[(*tri)[3 * triIndex + 1]].imag() - zc.imag()); 1021 DDouble y23 = (y2 - y3); 1022 DDouble x32 = (x3 - x2); 1023 DDouble y31 = (y3 - y1); 1024 DDouble x13 = (x1 - x3); 1025 DDouble det = (y23 * x13 - x32 * y31); 1026 DDouble minD = (det<0)?det:0 ; // we MUST be agnostic as to the orientation of the triangle, even if TRIANGULATE is clockwise. 1027 DDouble maxD = (det>0)?det:0; 1028 for (SizeT iy = pymin; iy <= pymax; ++iy) { 1029 DDouble dy = y[iy] - y3; 1030 DLong start = iy*nx; 1031 for (SizeT ix = pxmin; ix <= pxmax; ++ix) { 1032 DDouble dx = x[ix] - x3; 1033 DDouble a = y23 * dx + x32 * dy; 1034 if (a < minD || a > maxD) continue; 1035 DDouble b = y31 * dx + x13 * dy; 1036 if (b < minD || b > maxD) continue; 1037 DDouble c = det - a - b; 1038 if (c < minD || c > maxD) continue; 1039 DDouble valbisR = (a * zac.real() + b * zbc.real()) / det + zc.real(); 1040 DDouble valbisI = (a * zac.imag() + b * zbc.imag()) / det + zc.imag(); 1041 if ((dominvalue && valbisR < minVal.real()) || (domaxvalue && valbisR > maxVal.real())) { 1042 valbisR = (*res)[start + ix].real(); 1043 } 1044 if ((dominvalue && valbisI < minVal.imag()) || (domaxvalue && valbisI > maxVal.imag())) { 1045 valbisI = (*res)[start + ix].imag(); 1046 } 1047 (*res)[start + ix] = std::complex<double>(valbisR, valbisI); 1048 } //ix loop 1049 } // iy 1050 } //tri 1051 } else { 1052 for (SizeT triIndex = 0; triIndex < ntri; ++triIndex) { 1053 // first, find whose region of grid is concerned by this triangle: boundingBox. 1054 DDouble xmin; 1055 DDouble xmax; 1056 DDouble ymin; 1057 DDouble ymax; 1058 DDouble x1, y1, x2, y2, x3, y3; 1059 x1 = (*xVal)[(*tri)[3 * triIndex]]; 1060 y1 = (*yVal)[(*tri)[3 * triIndex]]; 1061 x2 = (*xVal)[(*tri)[3 * triIndex + 1]]; 1062 y2 = (*yVal)[(*tri)[3 * triIndex + 1]]; 1063 x3 = (*xVal)[(*tri)[3 * triIndex + 2]]; 1064 y3 = (*yVal)[(*tri)[3 * triIndex + 2]]; 1065 xmin = x1 < x2 ? (x1 < x3 ? x1 : x3) : (x2 < x3 ? x2 : x3); 1066 ymin = y1 < y2 ? (y1 < y3 ? y1 : y3) : (y2 < y3 ? y2 : y3); 1067 xmax = x1 > x2 ? (x1 > x3 ? x1 : x3) : (x2 > x3 ? x2 : x3); 1068 ymax = y1 > y2 ? (y1 > y3 ? y1 : y3) : (y2 > y3 ? y2 : y3); 1069 1070 // find first and last index of x and y that will intersect this BBox 1071 DLong k = 0; 1072 while (k < (nx - 1) && x[k + 1] < xmin) { 1073 k++; 1074 } 1075 int pxmin = k; 1076 while (k < (nx - 1) && x[k + 1] < xmax) { 1077 k++; 1078 } 1079 int pxmax = k; 1080 k = 0; 1081 while (k < (ny - 1) && y[k + 1] < ymin) { 1082 k++; 1083 } 1084 int pymin = k; 1085 while (k < (ny - 1) && y[k + 1] < ymax) { 1086 k++; 1087 } 1088 int pymax = k; 1089 1090 DComplexDbl zc = (*zVal)[(*tri)[3 * triIndex + 2]]; 1091 DComplexDbl zac = std::complex<double>((*zVal)[(*tri)[3 * triIndex]].real() - zc.real(), (*zVal)[(*tri)[3 * triIndex]].imag() - zc.imag()); 1092 DComplexDbl zbc = std::complex<double>((*zVal)[(*tri)[3 * triIndex + 1]].real() - zc.real(), (*zVal)[(*tri)[3 * triIndex + 1]].imag() - zc.imag()); 1093 DDouble y23 = (y2 - y3); 1094 DDouble x32 = (x3 - x2); 1095 DDouble y31 = (y3 - y1); 1096 DDouble x13 = (x1 - x3); 1097 DDouble det = (y23 * x13 - x32 * y31); 1098 DDouble minD = (det<0)?det:0 ; // we MUST be agnostic as to the orientation of the triangle, even if TRIANGULATE is clockwise. 1099 DDouble maxD = (det>0)?det:0; 1100 for (SizeT iy = pymin; iy <= pymax; ++iy) { 1101 DDouble dy = y[iy] - y3; 1102 DLong start = iy*nx; 1103 for (SizeT ix = pxmin; ix <= pxmax; ++ix) { 1104 DDouble dx = x[ix] - x3; 1105 DDouble a = y23 * dx + x32 * dy; 1106 if (a < minD || a > maxD) continue; 1107 DDouble b = y31 * dx + x13 * dy; 1108 if (b < minD || b > maxD) continue; 1109 DDouble c = det - a - b; 1110 if (c < minD || c > maxD) continue; 1111 DDouble valbisR = (a * zac.real() + b * zbc.real()) / det + zc.real(); 1112 DDouble valbisI = (a * zac.imag() + b * zbc.imag()) / det + zc.imag(); 1113 (*res)[start + ix] = std::complex<double>(valbisR, valbisI); 1114 } //ix loop 1115 } // iy 1116 } //tri 1117 } 1118 } 1119 trigrid_fun_plane(EnvT * e)1120 BaseGDL* trigrid_fun_plane(EnvT* e) { 1121 SizeT nParam = e->NParam(); 1122 if (nParam < 4) e->Throw("Incorrect number of arguments."); 1123 //OK, trigrid does not care if more than 3 args in sphere mode. Limit at 6 is done at interpreter level. 1124 1125 1126 BaseGDL* p0 = e->GetParDefined(0); //x 1127 BaseGDL* p1 = e->GetParDefined(1); //y 1128 BaseGDL* p2 = e->GetParDefined(2); //z 1129 1130 DType type=GDL_FLOAT; 1131 bool isComplex=(p2->Type()==GDL_COMPLEX || p2->Type()==GDL_COMPLEXDBL); 1132 //Double output is not only defined by z but also by x and y. 1133 if (p2->Type()==GDL_DOUBLE || p2->Type()==GDL_COMPLEXDBL) type=GDL_DOUBLE; 1134 else if (p1->Type()==GDL_DOUBLE || p1->Type()==GDL_COMPLEXDBL) type=GDL_DOUBLE; 1135 else if (p0->Type()==GDL_DOUBLE || p0->Type()==GDL_COMPLEXDBL) type=GDL_DOUBLE; 1136 1137 BaseGDL* p3 = e->GetParDefined(3); //triangles 1138 1139 if (p0->N_Elements() != p1->N_Elements() || 1140 p0->N_Elements() != p2->N_Elements() || 1141 p1->N_Elements() != p2->N_Elements()) 1142 e->Throw("X, Y, or Z array dimensions are incompatible."); 1143 1144 if (p0->N_Elements() < 3) e->Throw("Value of Bounds is out of allowed range."); 1145 1146 if (p3->Rank() == 0) e->Throw("Expression must be an array in this context: " + e->GetParString(0)); 1147 if (p3->N_Elements() % 3 != 0) e->Throw("Array of triangles incorrectly dimensioned."); 1148 1149 DLong ntri = p3->N_Elements() / 3; 1150 1151 if (p0->Rank() == 0) e->Throw("Expression must be an array in this context: " + e->GetParString(0)); 1152 if (p1->Rank() == 0) e->Throw("Expression must be an array in this context: " + e->GetParString(1)); 1153 if (p2->Rank() == 0) e->Throw("Expression must be an array in this context: " + e->GetParString(2)); 1154 1155 static int xoutIx=e->KeywordIx( "XOUT"); 1156 bool doXout=(e->KeywordSet(xoutIx)); 1157 static int youtIx=e->KeywordIx( "YOUT"); 1158 bool doYout=(e->KeywordSet(youtIx)); 1159 if ((doXout && !doYout) || (doYout && !doXout)) e->Throw("Incorrect number of arguments."); 1160 static int xgridIx=e->KeywordIx( "XGRID"); 1161 bool doXgrid=(e->KeywordPresent(xgridIx)); 1162 static int ygridIx=e->KeywordIx( "YGRID"); 1163 bool doYgrid=(e->KeywordPresent(ygridIx)); 1164 1165 static int inputIx=e->KeywordIx( "INPUT"); 1166 bool doInput=(e->KeywordPresent(inputIx)); 1167 1168 // Get NX, NY values if present 1169 DLong nx = 51; 1170 DLong ny = 51; 1171 static int nxIx = e->KeywordIx("NX"); 1172 static int nyIx = e->KeywordIx("NY"); 1173 bool canUseLimitsx=true; 1174 bool canUseLimitsy=true; 1175 if (e->KeywordSet(nxIx)) { canUseLimitsx=false; e->AssureLongScalarKW(nxIx, nx);} 1176 if (e->KeywordSet(nyIx)) { canUseLimitsy=false; e->AssureLongScalarKW(nyIx, ny);} 1177 if (nx < 1 || ny < 1) e->Throw("Array dimensions must be greater than 0."); 1178 // we will further on define GS, whose value may overwrite the nx, ny. 1179 1180 DDoubleGDL* GS = NULL; 1181 DDoubleGDL* limits = NULL; 1182 if (nParam > 4) 1183 { 1184 BaseGDL* p4 = e->GetParDefined(4); 1185 if (p4->Rank() == 0) e->Throw("Expression must be an array in this context: " + e->GetParString(4)); 1186 if (p4->N_Elements() != 2) e->Throw("Array must have 2 elements: " + e->GetParString(4)); 1187 GS = static_cast<DDoubleGDL*>(p4->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 1188 } 1189 if (nParam == 6) 1190 { 1191 BaseGDL* p5 = e->GetParDefined(5); 1192 if (p5->Rank() == 0) e->Throw("Expression must be an array in this context: " + e->GetParString(5)); 1193 if (p5->N_Elements() != 4) e->Throw("Array must have 4 elements: " + e->GetParString(5)); 1194 limits = static_cast<DDoubleGDL*> (p5->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 1195 } 1196 DLong minxEl; 1197 DLong maxxEl; 1198 DLong minyEl; 1199 DLong maxyEl; 1200 DDoubleGDL* xVal = static_cast<DDoubleGDL*>(p0->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 1201 DDoubleGDL* yVal = static_cast<DDoubleGDL*>(p1->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 1202 DLongGDL* tri = static_cast<DLongGDL*>(p3->Convert2(GDL_LONG, BaseGDL::COPY)); 1203 //get min max X Y. 1204 xVal->MinMax(&minxEl, &maxxEl, NULL, NULL, true); 1205 yVal->MinMax(&minyEl, &maxyEl, NULL, NULL, true); 1206 // Determine grid range 1207 DLong xref=0, yref=0; 1208 DDouble xval = (*xVal)[minxEl]; 1209 DDouble yval = (*yVal)[minyEl]; 1210 DDouble xrange = (*xVal)[maxxEl] - (*xVal)[minxEl]; 1211 DDouble yrange = (*yVal)[maxyEl] - (*yVal)[minyEl]; 1212 //compute World positions of each pixels. 1213 if (limits != NULL) 1214 { 1215 if (canUseLimitsx) { 1216 xval = (*limits)[0]; 1217 xrange = (*limits)[2] - xval; 1218 } 1219 if (canUseLimitsy) { 1220 yval = (*limits)[1]; 1221 yrange = (*limits)[3] - yval; 1222 } 1223 } 1224 // Determine grid spacing 1225 DDouble xinc = xrange / (nx-1); 1226 DDouble yinc = yrange / (ny-1); 1227 if (GS != NULL && canUseLimitsx) 1228 { 1229 xinc = (*GS)[0]; 1230 nx = (DLong) ceil(xrange / xinc) +1; 1231 } 1232 if (GS != NULL && canUseLimitsy) 1233 { 1234 yinc = (*GS)[1]; 1235 ny = (DLong) ceil(yrange / yinc) +1; 1236 } 1237 1238 DDouble *x, *y; 1239 bool cleanx=false; 1240 bool cleany=false; 1241 1242 if (doXout) { 1243 DDoubleGDL* xOut=e->GetKWAs<DDoubleGDL>(xoutIx); 1244 nx=xOut->N_Elements(); 1245 if (nx < 1) e->Throw("XOUT must have >1 elements."); 1246 x=(DDouble*)xOut->DataAddr(); 1247 } else { 1248 x = (DDouble*)malloc(nx*sizeof(DDouble)); 1249 for (int i=0; i<nx; ++i) x[i]=(i-xref)*xinc+xval; 1250 cleanx=true; 1251 } 1252 1253 if (doYout) { 1254 DDoubleGDL* yOut=e->GetKWAs<DDoubleGDL>(youtIx); 1255 ny=yOut->N_Elements(); 1256 if (ny < 1) e->Throw("YOUT must have >1 elements."); 1257 y=(DDouble*)yOut->DataAddr(); 1258 } else { 1259 y = (DDouble*)malloc(ny*sizeof(DDouble)); 1260 for (int j=0; j<ny; ++j) y[j]=(j-yref)*yinc+yval; 1261 cleany=true; 1262 } 1263 1264 //must be done here 1265 if (doXgrid) { 1266 DDoubleGDL *xgrid=new DDoubleGDL(nx, BaseGDL::NOZERO); 1267 for (int i=0; i<nx; ++i) (*xgrid)[i]=x[i]; 1268 e->SetKW(xgridIx, xgrid); 1269 } 1270 if (doYgrid) { 1271 DDoubleGDL *ygrid=new DDoubleGDL(ny, BaseGDL::NOZERO); 1272 for (int i=0; i<ny; ++i) (*ygrid)[i]=y[i]; 1273 e->SetKW(ygridIx, ygrid); 1274 } 1275 1276 // Setup return array 1277 DLong dims[2]; 1278 dims[0] = nx; 1279 dims[1] = ny; 1280 dimension dim((DLong *) dims, 2); 1281 1282 // if INPUT kw, test if input is an existing variable 1283 BaseGDL* inputKwData=NULL; 1284 DDoubleGDL* inputArray=NULL; 1285 bool inputArrayPresent=false; 1286 if (doInput) { 1287 inputKwData=e->GetKW(inputIx); 1288 if (!(inputKwData==NULL)) { //does exist, check compatibility 1289 if (inputKwData->Dim() != dim) e->Throw("Array has a corrupted descriptor"); 1290 inputArrayPresent=true; 1291 } 1292 } 1293 1294 if (isComplex) { 1295 DComplexDblGDL* minValG=NULL; 1296 static int minvalueIx=e->KeywordIx( "MIN_VALUE"); 1297 bool dominvalue=(e->KeywordPresent(minvalueIx)); 1298 if (dominvalue) minValG=e->GetKWAs<DComplexDblGDL>(minvalueIx); 1299 DComplexDblGDL* maxValG=NULL; 1300 static int maxvalueIx=e->KeywordIx( "MAX_VALUE"); 1301 bool domaxvalue=(e->KeywordPresent(maxvalueIx)); 1302 if (domaxvalue) maxValG=e->GetKWAs<DComplexDblGDL>(maxvalueIx); 1303 DComplexDblGDL* missValG=NULL; 1304 static int missvalueIx=e->KeywordIx( "MISSING"); 1305 bool doMiss=e->KeywordPresent(missvalueIx); 1306 if (doMiss) missValG=e->GetKWAs<DComplexDblGDL>(missvalueIx); 1307 DComplexDblGDL* zVal = static_cast<DComplexDblGDL*>(p2->Convert2(GDL_COMPLEXDBL, BaseGDL::COPY)); 1308 DComplexDblGDL* res; 1309 if (inputArrayPresent) res=static_cast<DComplexDblGDL*>(inputKwData->Convert2(GDL_COMPLEXDBL, BaseGDL::COPY)); 1310 else if (doMiss) { 1311 res = missValG->New(dim, BaseGDL::INIT); 1312 } else { 1313 res= new DComplexDblGDL(dim, BaseGDL::ZERO); 1314 } 1315 1316 DComplexDbl minVal=std::complex<double>(0,0); 1317 if (minValG!=NULL) minVal=std::complex<double>((*minValG)[0].real(),(*minValG)[0].imag()); 1318 DComplexDbl maxVal=std::complex<double>(0,0); 1319 if (maxValG!=NULL) maxVal=std::complex<double>((*maxValG)[0].real(),(*maxValG)[0].imag()); 1320 if (doXout) gdlGrid2DDataCpx(nx, x, ny, y, ntri, tri, xVal, yVal, zVal, domaxvalue, dominvalue, maxVal, minVal, res, inputArrayPresent); 1321 else gdlGrid2DDataCpx(nx, x, xref,xval,xinc, ny, y, yref,yval,yinc,ntri, tri, xVal, yVal, zVal, domaxvalue, dominvalue, maxVal, minVal, res, inputArrayPresent); 1322 if (type==GDL_FLOAT) return res->Convert2(GDL_COMPLEX,BaseGDL::CONVERT); else return res; 1323 }else{ 1324 DDouble minVal=0; 1325 static int minvalueIx=e->KeywordIx( "MIN_VALUE"); 1326 bool dominvalue=(e->KeywordPresent(minvalueIx)); 1327 if (dominvalue) e->AssureDoubleScalarKW(minvalueIx, minVal); 1328 DDouble maxVal=0; 1329 static int maxvalueIx=e->KeywordIx( "MAX_VALUE"); 1330 bool domaxvalue=(e->KeywordPresent(maxvalueIx)); 1331 if (domaxvalue) e->AssureDoubleScalarKW(maxvalueIx, maxVal); 1332 DDoubleGDL* missVal=NULL; 1333 static int missvalueIx=e->KeywordIx( "MISSING"); 1334 bool doMiss=e->KeywordPresent(missvalueIx); 1335 if (doMiss) missVal=e->GetKWAs<DDoubleGDL>(missvalueIx); 1336 DDoubleGDL* zVal = static_cast<DDoubleGDL*>(p2->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 1337 DDoubleGDL* res; 1338 if (inputArrayPresent) { 1339 res = static_cast<DDoubleGDL*> (inputKwData->Convert2(GDL_DOUBLE, BaseGDL::COPY)); 1340 } else if (doMiss) { 1341 res = missVal->New(dim, BaseGDL::INIT); 1342 } else { 1343 res = new DDoubleGDL(dim, BaseGDL::ZERO); 1344 } 1345 if (doXout) gdlGrid2DData< DDoubleGDL, DDouble>(nx, x, ny, y, ntri, tri, xVal, yVal, zVal, domaxvalue, dominvalue, maxVal, minVal, res, inputArrayPresent); 1346 else gdlGrid2DData< DDoubleGDL, DDouble>(nx, x, xref,xval,xinc, ny, y, yref,yval,yinc,ntri, tri, xVal, yVal, zVal, domaxvalue, dominvalue, maxVal, minVal, res, inputArrayPresent); 1347 if (type==GDL_FLOAT) return res->Convert2(GDL_FLOAT,BaseGDL::CONVERT); else return res; 1348 } 1349 } 1350 1351 trigrid_fun(EnvT * e)1352 BaseGDL* trigrid_fun(EnvT* e) { 1353 //just send to normal or spherical based on presence of SPHERE argument 1354 static int sphereIx = e->KeywordIx( "SPHERE"); 1355 if(e->KeywordPresent(sphereIx)) return trigrid_fun_spherical(e); else return trigrid_fun_plane(e); 1356 } 1357 grid_input(EnvT * e)1358 void grid_input(EnvT* e) { 1359 e->Throw("Writing in progress."); 1360 } 1361 1362 //to be written and do not forget to uncomment QHULL in CMakeLists and config.h.cmake 1363 // see http://www.geom.umn.edu/software/qhull/. Used also with plplot. 1364 #ifdef HAVE_QHULL qhull(EnvT * e)1365 void qhull ( EnvT* e) 1366 { 1367 e->Throw("Please Write this function in GDL."); 1368 } 1369 1370 qgrid3_fun(EnvT * e)1371 BaseGDL* qgrid3_fun ( EnvT* e) 1372 { 1373 e->Throw("Please Write this function in GDL."); 1374 return NULL; 1375 } 1376 #endif 1377 } 1378 1379