1 /***************************************************************************
2 plotting.cpp - GDL routines for plotting
3 -------------------
4 begin : July 22 2002
5 copyright : (C) 2002-2011 by Marc Schellens et al.
6 email : m_schellens@users.sf.net
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 "dinterpreter.hpp"
20 #include "plotting.hpp"
21
22 #ifdef _MSC_VER
23 #define snprintf _snprintf
24 #endif
25
26 namespace lib
27 {
28
29 using namespace std;
30 // using std::isinf;
31 #ifdef _MSC_VER
32 #define finite _finite
33 #else
34 using std::isnan;
35 #endif
36
37 //static values
38 static DLong savedStyle=0;
39 static DDouble savedPointX=0.0;
40 static DDouble savedPointY=0.0;
41 static gdlSavebox saveBox;
42 static DFloat sym1x[6]={1, -1, 0, 0, 0, 0}; // +
43 static DFloat sym1y[6]={0, 0, 0, -1, 1, 0}; // +
44 static DFloat sym2x[12]= {1, -1, 0, 0, 0, 0,1,-1,0,1,-1, 0}; //*
45 static DFloat sym2y[12]= {0, 0, 0, -1, 1,0,1,-1,0,-1,1, 0}; // *
46 static DFloat sym3x[2]={0,0}; // dot. On PostScript device, x1=x1 and y2=1 creates a round dot.
47 static DFloat sym3y[2]={0,0}; // .
48 static DFloat sym4x[5]={ 0, 1, 0, -1, 0 }; //diamond.
49 static DFloat sym4y[5]={ 1, 0, -1, 0, 1 }; //diamond.
50 static DFloat sym5x[4]={ -1, 0, 1, -1 }; // triangle up.
51 static DFloat sym5y[4]={ -1, 1, -1, -1 }; // triangle up.
52 static DFloat sym6x[5]={ -1, 1, 1, -1, -1 }; //square
53 static DFloat sym6y[5]={ 1, 1, -1, -1, 1 }; //square
54 static DFloat sym7x[6]= {1,-1,0,1,-1,0}; //x
55 static DFloat sym7y[6]= {1,-1,0,-1,1,0}; //x
56 DLong syml[7]={6,12,2,5,4,5,6};
57 static DFloat sym3xalt[6]={-0.2,-0.2,0.2,0.2,-0.2,0}; // big dot.
58 static DFloat sym3yalt[6]={-0.2,0.2,0.2,-0.2,-0.2,0}; // .
59 DLong syml_alt=6;
60
61 struct LOCALUSYM {
62 DLong nusym;
63 DInt fill;
64 DFloat usymx[49];
65 DFloat usymy[49];
66 bool hasColor;
67 DLong color;
68 bool hasThick;
69 DFloat thick;
70 };
71 static LOCALUSYM localusym;
72
getSaveBox()73 gdlSavebox* getSaveBox(){return &saveBox;}
74
gdlDoRangeExtrema(DDoubleGDL * xVal,DDoubleGDL * yVal,DDouble & min,DDouble & max,DDouble xmin,DDouble xmax,bool doMinMax,DDouble minVal,DDouble maxVal)75 void gdlDoRangeExtrema(DDoubleGDL *xVal, DDoubleGDL *yVal, DDouble &min, DDouble &max, DDouble xmin, DDouble xmax, bool doMinMax, DDouble minVal, DDouble maxVal)
76 {
77 DDouble valx, valy;
78 SizeT i,k;
79 DLong n=xVal->N_Elements();
80 if(n!=yVal->N_Elements()) return;
81 for (i=0,k=0 ; i<n ; ++i)
82 {
83 //look only in range x=[xmin,xmax]
84 valx=(*xVal)[i];
85 if ( (valx<xmin || valx>xmax || isnan(valx))) continue;
86 else {
87 //min and max of y if not NaN and in range [minVal, maxVal] if doMinMax=yes (min_value, max_value keywords)
88 valy=(*yVal)[i];
89 if ((doMinMax && (valy<minVal || valy>maxVal )) || isnan(valy)) continue;
90 else {if(k==0) {min=valy; max=valy;} else {min=gdlPlot_Min(min,valy); max=gdlPlot_Max(max,valy);}}
91 }
92 k++;
93 }
94 }
95
GetMinMaxVal(DDoubleGDL * val,double * minVal,double * maxVal)96 void GetMinMaxVal(DDoubleGDL* val, double* minVal, double* maxVal)
97 {
98 #define UNDEF_RANGE_VALUE 1E-12
99 DLong minE, maxE;
100 const bool omitNaN=true;
101 val->MinMax(&minE, &maxE, NULL, NULL, omitNaN);
102 if ( minVal!=NULL ) {
103 *minVal=(*val)[ minE];
104 if (isnan(*minVal)) *minVal = UNDEF_RANGE_VALUE;
105 }
106 if ( maxVal!=NULL ) {
107 *maxVal=(*val)[ maxE];
108 if (isnan(*maxVal)) *maxVal = 1.0;
109 }
110 if ((*maxVal)==(*minVal)) *maxVal=*minVal+1.0;
111 #undef UNDEF_RANGE_VALUE
112 }
113
GetMinMaxValuesForSubset(DDoubleGDL * val,DDouble & minVal,DDouble & maxVal,SizeT FinalElement)114 void GetMinMaxValuesForSubset(DDoubleGDL* val, DDouble &minVal, DDouble &maxVal, SizeT FinalElement)
115 {
116 #define UNDEF_RANGE_VALUE 1E-12
117 DLong minE, maxE;
118 const bool omitNaN=true;
119 val->MinMax(&minE, &maxE, NULL, NULL, omitNaN, 0, FinalElement);
120 minVal=(*val)[ minE];
121 if (isnan(minVal)) minVal = UNDEF_RANGE_VALUE;
122 maxVal=(*val)[ maxE];
123 if (isnan(maxVal)) maxVal = 1.0;
124 if (maxVal==minVal) maxVal=minVal+1.0;
125 #undef UNDEF_RANGE_VALUE
126 }
127
AutoTick(DDouble x)128 PLFLT AutoTick(DDouble x)
129 {
130 if ( x==0.0 ) return 1.0;
131
132 DLong n=static_cast<DLong>(floor(log10(x/3.5)));
133 DDouble y=(x/(3.5*pow(10., static_cast<double>(n))));
134 DLong m=0;
135 if ( y>=1&&y<2 )
136 m=1;
137 else if ( y>=2&&y<5 )
138 m=2;
139 else if ( y>=5 )
140 m=5;
141
142 PLFLT intv=(PLFLT)(m*pow(10., static_cast<double>(n)));
143 return intv;
144 }
145
AutoIntv(DDouble x)146 PLFLT AutoIntv(DDouble x)
147 {
148 if ( x==0.0 )
149 {
150 // cout << "zero"<<endl;
151 return 1.0;
152 }
153
154 DLong n=static_cast<DLong>(floor(log10(x/2.82)));
155 DDouble y=(x/(2.82*pow(10., static_cast<double>(n))));
156 DLong m=0;
157 if ( y>=1&&y<2 )
158 m=1;
159 else if ( y>=2&&y<4.47 )
160 m=2;
161 else if ( y>=4.47 )
162 m=5;
163
164 // cout << "AutoIntv" << x << " " << y << endl;
165
166 PLFLT intv=(PLFLT)(m*pow(10., static_cast<double>(n)));
167 return intv;
168 }
169
170
171
172 #define EXTENDED_DEFAULT_LOGRANGE 12
173 //protect from (inverted, strange) axis log values
gdlHandleUnwantedAxisValue(DDouble & min,DDouble & max,bool log)174 void gdlHandleUnwantedAxisValue(DDouble &min, DDouble &max, bool log)
175 {
176 bool invert=FALSE;
177 DDouble val_min, val_max;
178 if (!log) return;
179
180 if(max-min >= 0)
181 {
182 val_min=min;
183 val_max=max;
184 invert=FALSE;
185 } else {
186 val_min=max;
187 val_max=min;
188 invert=TRUE;
189 }
190
191 if ( val_min<=0. )
192 {
193 if ( val_max<=0. )
194 {
195 val_min=-EXTENDED_DEFAULT_LOGRANGE;
196 val_max=0.0;
197 }
198 else
199 {
200 val_min=log10(val_max)-EXTENDED_DEFAULT_LOGRANGE;
201 val_max=log10(val_max);
202 }
203 }
204 else
205 {
206 val_min=log10(val_min);
207 val_max=log10(val_max);
208 }
209 if (invert)
210 {
211 min=pow(10.,val_max);
212 max=pow(10.,val_min);
213 } else {
214 min=pow(10.,val_min);
215 max=pow(10.,val_max);
216 }
217
218 }
219 #undef EXTENDED_DEFAULT_LOGRANGE
220
221 //improved version of "AutoIntv" for:
222 // 1/ better managing ranges when all the data have same value
223 // 2/ mimic IDL behavior when data are all positive
224 // please notice that (val_min, val_max) will be changed
225 // and "epsilon" is a coefficient if "extended range" is expected
226 // input: linear min and max, output: linear min and max.
227
228 // NOTE GD: this function should be rewritten, documented and tested correctly. Most often, the
229 // plots are not exactly what IDL does in the same conditions. The reasons for the choices should be
230 // clearly described in the code, to be checked by others.
231
gdlAdjustAxisRange(EnvT * e,int axisId,DDouble & start,DDouble & end,bool log,int code)232 PLFLT gdlAdjustAxisRange(EnvT* e, int axisId, DDouble &start, DDouble &end, bool log /* = false */, int code /* = 0 */) {
233 gdlHandleUnwantedAxisValue(start, end, log);
234
235 DDouble min, max;
236 bool invert = FALSE;
237
238 if (end - start >= 0) {
239 min = start;
240 max = end;
241 invert = FALSE;
242 } else { //never happens by construction!!!!!!!!!!!!!!!!!!!!!!!
243 min = end;
244 max = start;
245 invert = TRUE;
246 }
247
248 PLFLT intv = 1.;
249 int cas = 0;
250 DDouble x;
251 bool debug = false;
252 if (debug) {
253 cout << "init: " << min << " " << max << endl;
254 }
255 // case "all below ABS((MACHAR()).xmin)"
256 //needs example.
257
258 // if (!log && (abs(max) <= -std::numeric_limits<DDouble>::max())) {
259 // min = DDouble(0.);
260 // max = DDouble(1.);
261 // intv = (PLFLT) (2.);
262 // cas = 1;
263 // }
264
265 if (log) {
266 min = log10(min);
267 max = log10(max);
268 }
269
270 // case "all values are equal"
271 if (cas == 0) {
272 x = max - min;
273 if (abs(x) <= std::numeric_limits<DDouble>::min()) {
274 DDouble val_ref;
275 val_ref = max;
276 if (0.98 * min < val_ref) { // positive case
277 max = 1.02 * val_ref;
278 min = 0.98 * val_ref;
279 } else { // negative case
280 max = 0.98 * val_ref;
281 min = 1.02 * val_ref;
282 }
283 if (debug) {
284 cout << "Rescale : " << min << " " << max << endl;
285 }
286 }
287 }
288
289 // general case (only negative OR negative and positive)
290 if (cas == 0) //rounding is not aka idl due to use of ceil and floor. TBD.
291 {
292 x = max - min;
293 //correct this for calendar values (round to nearest year, month, etc)
294 if ( code > 0) {
295 if (code ==7 ) {
296 if(x>=366) code=1;
297 else if(x>=32) code=2;
298 else if(x>=1.1) code=3;
299 else if(x*24>=1.1) code=4;
300 else if(x*24*60>=1.1) code=5;
301 else code=6;
302 }
303 static int monthSize[]={31,28,31,30,31,30,31,31,30,31,30,31};
304 DLong Day1,Day2 , Year1,Year2 , Hour1,Hour2 , Minute1,Minute2, MonthNum1,MonthNum2;
305 DLong idow1,icap1,idow2,icap2;
306 DDouble Seconde1,Seconde2;
307 j2ymdhms(min, MonthNum1 , Day1 , Year1 , Hour1 , Minute1, Seconde1, idow1, icap1);
308 MonthNum1++; //j2ymdhms gives back Month number in the [0-11] range for indexing month name tables. pity.
309 j2ymdhms(max, MonthNum2 , Day2 , Year2 , Hour2 , Minute2, Seconde2, idow2, icap2);
310 MonthNum2++;
311 switch(code){
312 case 1:
313 // day mon year h m s.s
314 dateToJD(min, 1, 1, Year1, 0, 0, 0.0);
315 dateToJD(max, 1, 1, Year2+1, 0, 0, 0.0);
316 break;
317 case 2:
318 dateToJD(min, 1, MonthNum1, Year1, 0, 0, 0.0);
319 MonthNum2++;
320 if (MonthNum2 > 12) {MonthNum2-=12; Year2+=1;}
321 dateToJD(max, 1, MonthNum2, Year2, 0, 0, 0.0);
322 break;
323 case 3:
324 dateToJD(min, Day1, MonthNum1, Year1, 0, 0, 0.0);
325 Day2++;
326 if (Day2 > monthSize[MonthNum2]) {Day2-=monthSize[MonthNum2]; MonthNum2+=1;}
327 if (MonthNum2 > 12) {MonthNum2-=12; Year2+=1;}
328 dateToJD(max, Day2, MonthNum2, Year2, 0, 0, 0.0);
329 break;
330 case 4:
331 dateToJD(min, Day1, MonthNum1, Year1, Hour1, 0, 0.0);
332 Hour2++;
333 if (Hour2 > 23) {Hour2-=24; Day2+=1;}
334 if (Day2 > monthSize[MonthNum2]) {Day2-=monthSize[MonthNum2]; MonthNum2+=1;}
335 if (MonthNum2 > 12) {MonthNum2-=12; Year2+=1;}
336 dateToJD(max, Day2, MonthNum2, Year2, Hour2, 0, 0.0);
337 break;
338 case 5:
339 dateToJD(min, Day1, MonthNum1, Year1, Hour1, Minute1, 0.0);
340 Minute2++;
341 if (Minute2 > 59) {Minute2-=60; Hour2+=1;}
342 if (Hour2 > 23) {Hour2-=24; Day2+=1;}
343 if (Day2 > monthSize[MonthNum2]) {Day2-=monthSize[MonthNum2]; MonthNum2+=1;}
344 if (MonthNum2 > 12) {MonthNum2-=12; Year2+=1;}
345 dateToJD(max, Day2, MonthNum2, Year2, Hour2, Minute2, 0.0);
346 break;
347 case 6:
348 dateToJD(min, Day1, MonthNum1, Year1, Hour1, Minute1, Seconde1);
349 Seconde2++;
350 if (Seconde2 > 59) {Seconde2-=60; Minute2+=1;}
351 if (Minute2 > 59) {Minute2-=60; Hour2+=1;}
352 if (Hour2 > 23) {Hour2-=24; Day2+=1;}
353 if (Day2 > monthSize[MonthNum2]) {Day2-=monthSize[MonthNum2]; MonthNum2+=1;}
354 if (MonthNum2 > 12) {MonthNum2-=12; Year2+=1;}
355 dateToJD(max, Day2, MonthNum2, Year2, Hour2, Minute2, Seconde2);
356 break;
357 default:
358 break;
359 }
360 }
361 else {
362 intv = AutoIntv(x);
363 if (log) {
364 max = ceil((max / intv) * intv);
365 min = floor((min / intv) * intv);
366 } else {
367 max = ceil(max / intv) * intv;
368 min = floor(min / intv) * intv;
369 }
370 }
371
372
373 }
374
375 if (debug) {
376 cout << "cas: " << cas << " new range: " << min << " " << max << endl;
377 }
378 //give back non-log values
379 if (log) {
380 min = pow(10, min);
381 max = pow(10, max);
382 }
383
384 //check if tickinterval would make more than 59 ticks (IDL apparent limit). In which case, IDL plots only the first 59 intervals:
385 DDouble TickInterval;
386 gdlGetDesiredAxisTickInterval(e, axisId, TickInterval);
387 if ( TickInterval > 0.0 ) if ((max-min)/TickInterval > 59) max=min+59.0*TickInterval;
388
389 if (invert) {
390 start = max;
391 end = min;
392 } else {
393 start = min;
394 end = max;
395 }
396
397
398 return intv;
399 }
400
CheckMargin(GDLGStream * actStream,DFloat xMarginL,DFloat xMarginR,DFloat yMarginB,DFloat yMarginT,PLFLT & xMR,PLFLT & xML,PLFLT & yMB,PLFLT & yMT)401 void CheckMargin(GDLGStream* actStream,
402 DFloat xMarginL,
403 DFloat xMarginR,
404 DFloat yMarginB,
405 DFloat yMarginT,
406 PLFLT& xMR,
407 PLFLT& xML,
408 PLFLT& yMB,
409 PLFLT& yMT)
410 {
411 PLFLT sclx=actStream->dCharLength()/actStream->xSubPageSize(); //current char length/subpage size
412 xML=xMarginL*sclx; //margin as percentage of subpage
413 xMR=xMarginR*sclx;
414 PLFLT scly=actStream->dLineSpacing()/actStream->ySubPageSize(); //current char length/subpage size
415 yMB=(yMarginB)*scly;
416 yMT=(yMarginT)*scly; //to allow subscripts and superscripts (as in IDL)
417
418 if ( xML+xMR>=1.0 )
419 {
420 Message("XMARGIN to large (adjusted).");
421 PLFLT xMMult=xML+xMR;
422 xML/=xMMult*1.5;
423 xMR/=xMMult*1.5;
424 }
425 if ( yMB+yMT>=1.0 )
426 {
427 Message("YMARGIN to large (adjusted).");
428 PLFLT yMMult=yMB+yMT;
429 yMB/=yMMult*1.5;
430 yMT/=yMMult*1.5;
431 }
432 }
433
setIsoPort(GDLGStream * actStream,PLFLT x1,PLFLT x2,PLFLT y1,PLFLT y2,PLFLT aspect)434 void setIsoPort(GDLGStream* actStream,
435 PLFLT x1,
436 PLFLT x2,
437 PLFLT y1,
438 PLFLT y2,
439 PLFLT aspect)
440 {
441 PLFLT X1, X2, Y1, Y2, X1s, X2s, Y1s, Y2s, displacx, displacy, scalex, scaley, offsetx, offsety;
442 if ( aspect<=0.0 )
443 {
444 actStream->vpor(x1, x2, y1, y2);
445 return;
446 }
447 // here we need too compensate for the change of aspect due to eventual !P.MULTI plots
448 actStream->vpor(x1, x2, y1, y2); //ask for non-iso window
449 actStream->gvpd(X1, X2, Y1, Y2); //get viewport values
450 //compute relation desiredViewport-page viewport x=scalex*X+offsetx:
451 scalex=(x2-x1)/(X2-X1);
452 offsetx=(x1*X2-x2*X1)/(X2-X1);
453 scaley=(y2-y1)/(Y2-Y1);
454 offsety=(y1*Y2-y2*Y1)/(Y2-Y1);
455 //ask for wiewport scaled to isotropic by plplot
456 actStream->vpas(x1, x2, y1, y2, aspect);
457 //retrieve values
458 actStream->gvpd(X1s, X2s, Y1s, Y2s);
459 //measure displacement
460 displacx=X1s-X1;
461 displacy=Y1s-Y1;
462 //set wiewport scaled by plplot, displaced, as vpor using above linear transformation
463 x1=(X1s-displacx)*scalex+offsetx;
464 x2=(X2s-displacx)*scalex+offsetx;
465 y1=(Y1s-displacy)*scaley+offsety;
466 y2=(Y2s-displacy)*scaley+offsety;
467 actStream->vpor(x1, x2, y1, y2);
468 }
469
470
471
GetSFromPlotStructs(DDouble ** sx,DDouble ** sy,DDouble ** sz)472 void GetSFromPlotStructs(DDouble **sx, DDouble **sy, DDouble **sz)
473 {
474 DStructGDL* xStruct=SysVar::X(); //MUST NOT BE STATIC, due to .reset
475 DStructGDL* yStruct=SysVar::Y(); //MUST NOT BE STATIC, due to .reset
476 DStructGDL* zStruct=SysVar::Z(); //MUST NOT BE STATIC, due to .reset
477 unsigned sxTag=xStruct->Desc()->TagIndex("S");
478 unsigned syTag=yStruct->Desc()->TagIndex("S");
479 unsigned szTag=zStruct->Desc()->TagIndex("S");
480 if (sx != NULL) *sx= &(*static_cast<DDoubleGDL*>(xStruct->GetTag(sxTag, 0)))[0];
481 if (sy != NULL) *sy= &(*static_cast<DDoubleGDL*>(yStruct->GetTag(syTag, 0)))[0];
482 if (sz != NULL) *sz= &(*static_cast<DDoubleGDL*>(zStruct->GetTag(szTag, 0)))[0];
483 }
484
GetWFromPlotStructs(DFloat ** wx,DFloat ** wy)485 void GetWFromPlotStructs(DFloat **wx, DFloat **wy)
486 {
487 DStructGDL* xStruct=SysVar::X(); //MUST NOT BE STATIC, due to .reset
488 DStructGDL* yStruct=SysVar::Y(); //MUST NOT BE STATIC, due to .reset
489 unsigned xwindowTag=xStruct->Desc()->TagIndex("WINDOW");
490 unsigned ywindowTag=yStruct->Desc()->TagIndex("WINDOW");
491 *wx= &(*static_cast<DFloatGDL*>(xStruct->GetTag(xwindowTag, 0)))[0];
492 *wy= &(*static_cast<DFloatGDL*>(yStruct->GetTag(ywindowTag, 0)))[0];
493 }
setPlplotScale(GDLGStream * a)494 void setPlplotScale(GDLGStream* a)
495 {
496 DDouble *sx, *sy;
497 GetSFromPlotStructs( &sx, &sy );
498
499 DDouble xStart, xEnd, yStart, yEnd;
500 xStart=-sx[0]/sx[1];
501 xEnd=(1-sx[0])/sx[1];
502 yStart=-sy[0]/sy[1];
503 yEnd=(1-sy[0])/sy[1];
504 a->wind(xStart, xEnd, yStart, yEnd);
505 }
DataCoordLimits(DDouble * sx,DDouble * sy,DFloat * wx,DFloat * wy,DDouble * xStart,DDouble * xEnd,DDouble * yStart,DDouble * yEnd,bool clip_by_default)506 void DataCoordLimits(DDouble *sx, DDouble *sy, DFloat *wx, DFloat *wy,
507 DDouble *xStart, DDouble *xEnd, DDouble *yStart, DDouble *yEnd, bool clip_by_default)
508 {
509 *xStart=(wx[0]-sx[0])/sx[1];
510 *xEnd=(wx[1]-sx[0])/sx[1];
511 *yStart=(wy[0]-sy[0])/sy[1];
512 *yEnd=(wy[1]-sy[0])/sy[1];
513
514 // patch from Joanna (tracker item no. 3029409, see test_clip.pro)
515 if ( !clip_by_default )
516 {
517 // cout << "joanna" << endl;
518 DFloat wxlen=wx[1]-wx[0];
519 DFloat wylen=wy[1]-wy[0];
520 DFloat xlen= *xEnd- *xStart;
521 DFloat ylen= *yEnd- *yStart;
522 *xStart= *xStart-xlen/wxlen*wx[0];
523 *xEnd= *xEnd+xlen/wxlen*(1-wx[1]);
524 *yStart= *yStart-ylen/wylen*wy[0];
525 *yEnd= *yEnd+ylen/wylen*(1-wy[1]);
526 }
527 // cout << *xStart <<" "<< *xEnd << " "<< *yStart <<" "<< *yEnd << ""<< endl;
528 }
529
GetUsym(DLong ** n,DInt ** do_fill,DFloat ** x,DFloat ** y,bool ** do_color,DLong ** usymColor,bool ** do_thick,DFloat ** usymThick)530 void GetUsym(DLong **n, DInt **do_fill, DFloat **x, DFloat **y, bool **do_color, DLong **usymColor , bool **do_thick, DFloat **usymThick)
531 {
532 *n= &(localusym.nusym);
533 *do_fill= &(localusym.fill);
534 *do_color= &(localusym.hasColor);
535 *do_thick= &(localusym.hasThick);
536 *usymColor=&(localusym.color);
537 *usymThick=&(localusym.thick);
538 *x=localusym.usymx;
539 *y=localusym.usymy;
540 }
541
SetUsym(DLong n,DInt do_fill,DFloat * x,DFloat * y,bool usersymhascolor,DLong usymColor,bool usersymhasthick,DFloat usymThick)542 void SetUsym(DLong n, DInt do_fill, DFloat *x, DFloat *y, bool usersymhascolor, DLong usymColor , bool usersymhasthick, DFloat usymThick )
543 {
544 localusym.nusym=n;
545 localusym.fill=do_fill;
546 for ( int i=0; i<n; i++ )
547 {
548 localusym.usymx[i]=x[i];
549 localusym.usymy[i]=y[i];
550 }
551 localusym.hasColor=usersymhascolor;
552 localusym.hasThick=usersymhasthick;
553 localusym.color=usymColor;
554 localusym.thick=usymThick;
555 }
556
557 //This is the good way to get world start end end values.
GetCurrentUserLimits(GDLGStream * a,DDouble & xStart,DDouble & xEnd,DDouble & yStart,DDouble & yEnd)558 void GetCurrentUserLimits(GDLGStream *a, DDouble &xStart, DDouble &xEnd, DDouble &yStart, DDouble &yEnd)
559 {
560 DDouble *sx, *sy;
561 GetSFromPlotStructs( &sx, &sy );
562 DFloat *wx, *wy;
563 GetWFromPlotStructs(&wx, &wy);
564 PLFLT x1,x2,y1,y2;
565 x1=wx[0];x2=wx[1];y1=wy[0];y2=wy[1];
566 xStart=(x1-sx[0])/sx[1];
567 xEnd=(x2-sx[0])/sx[1];
568 yStart=(y1-sy[0])/sy[1];
569 yEnd=(y2-sy[0])/sy[1];
570 //probably overkill now...
571 if ((yStart == yEnd) || (xStart == xEnd))
572 {
573 if (yStart != 0.0 && yStart == yEnd){
574 Message("PLOTS: !Y.CRANGE ERROR, setting to [0,1]");
575 yStart = 0;
576 yEnd = 1;
577 }
578 if (xStart != 0.0 && xStart == xEnd){
579 Message("PLOTS: !X.CRANGE ERROR, setting to [0,1]");
580 xStart = 0;
581 xEnd = 1;
582 }
583 }
584 }
585
ac_histo(GDLGStream * a,int i_buff,PLFLT * x_buff,PLFLT * y_buff,bool xLog)586 void ac_histo(GDLGStream *a, int i_buff, PLFLT *x_buff, PLFLT *y_buff, bool xLog)
587 {
588 PLFLT x, x1, y, y1, val;
589 for ( int jj=1; jj<i_buff; ++jj )
590 {
591 x1=x_buff[jj-1];
592 x=x_buff[jj];
593 y1=y_buff[jj-1];
594 y=y_buff[jj];
595 // cf patch 3567803
596 if ( xLog )
597 {
598 // val=log10((pow(10.0,x1)+pow(10.0,x))/2.0);
599 val=x1+log10(0.5+0.5*(pow(10.0, x-x1)));
600 }
601 else
602 {
603 val=(x1+x)/2.0;
604 }
605 a->join(x1, y1, val, y1);
606 a->join(val, y1, val, y);
607 a->join(val, y, x, y);
608 }
609 }
610
611
stopClipping(GDLGStream * a)612 void stopClipping(GDLGStream *a)
613 {
614 if (saveBox.initialized) {
615 a->vpor(saveBox.nx1, saveBox.nx2, saveBox.ny1, saveBox.ny2); //restore norm of current box
616 a->wind(saveBox.wx1, saveBox.wx2, saveBox.wy1, saveBox.wy2); //give back world of current box
617 } else cerr<<"plot \"savebox\" not initialized, please report" <<endl;
618 }
619
saveLastPoint(GDLGStream * a,DDouble wx,DDouble wy)620 void saveLastPoint(GDLGStream *a, DDouble wx, DDouble wy)
621 {
622 a->WorldToNormedDevice(wx, wy, savedPointX, savedPointY);
623 if (GDL_DEBUG_PLSTREAM) fprintf(stderr,"saveLastPoint as %lf %lf\n", savedPointX, savedPointY);
624 }
625
getLastPoint(GDLGStream * a,DDouble & wx,DDouble & wy)626 void getLastPoint(GDLGStream *a, DDouble& wx, DDouble& wy)
627 {
628 a->NormedDeviceToWorld(savedPointX, savedPointY, wx, wy);
629 if (GDL_DEBUG_PLSTREAM) fprintf(stderr,"getLastPoint: Got dev: %lf %lf giving %lf %lf world\n", savedPointX, savedPointY, wx, wy);
630 }
631 //LINESTYLE
gdlLineStyle(GDLGStream * a,DLong style)632 void gdlLineStyle(GDLGStream *a, DLong style)
633 {
634 //set saved Satle to nex style:
635 savedStyle=style;
636 static PLINT mark1[]={75};
637 static PLINT space1[]={1500};
638 static PLINT mark2[]={1500};
639 static PLINT space2[]={1500};
640 static PLINT mark3[]={1500, 100};
641 static PLINT space3[]={1000, 1000};
642 static PLINT mark4[]={1500, 100, 100, 100};
643 static PLINT space4[]={1000, 1000, 1000, 1000};
644 static PLINT mark5[]={3000};
645 static PLINT space5[]={1500}; // see plplot-5.5.3/examples/c++/x09.cc
646 switch(style)
647 {
648 case 0:
649 a->styl(0, mark1, space1);
650 return;
651 case 1:
652 a->styl(1, mark1, space1);
653 return;
654 case 2:
655 a->styl(1, mark2, space2);
656 return;
657 case 3:
658 a->styl(2, mark3, space3);
659 return;
660 case 4:
661 a->styl(4, mark4, space4);
662 return;
663 case 5:
664 a->styl(1, mark5, space5);
665 return;
666 default:
667 a->styl(0, NULL, NULL);
668 return;
669 }
670 }
gdlGetCurrentStyle()671 DLong gdlGetCurrentStyle(){
672 return savedStyle;
673 }
674
675 ///
676 /// Draws a line along xVal, yVal
677 /// @param general environnement pointer
678 /// @param graphic stream
679 /// @param xVal pointer on DDoubleGDL x values
680 /// @param yVal pointer on DDoubleGDL y values
681 /// @param minVal DDouble min value to plot.
682 /// @param maxVal DDouble max value to plot.
683 /// @param doMinMax bool do we use minval & maxval above?
684 /// @param xLog bool scale is log in x
685 /// @param yLog bool scale is log in y
686 /// @param psym DLong plotting symbol code
687 /// @param append bool values must be drawn starting from last plotted value
688 /// @param color DLongGDL* pointer to color list (NULL if no use)
689 ///
690
draw_polyline(GDLGStream * a,DDoubleGDL * xVal,DDoubleGDL * yVal,DDouble minVal,DDouble maxVal,bool doMinMax,bool xLog,bool yLog,DLong psym,bool useProjInfo,bool append,DLongGDL * color)691 void draw_polyline(GDLGStream *a,
692 DDoubleGDL *xVal, DDoubleGDL *yVal,
693 DDouble minVal, DDouble maxVal, bool doMinMax,
694 bool xLog, bool yLog, //end non-implicit parameters
695 DLong psym, bool useProjInfo, bool append, DLongGDL *color)
696 {
697 bool docolor=(color != NULL);
698
699 // Get decomposed value for colors
700 DLong decomposed=GraphicsDevice::GetDevice()->GetDecomposed();
701 //if docolor, do we really have more than one color?
702 if (docolor) if (color->N_Elements() == 1) { //do the job once and forget it after.
703 a->Color ( ( *color )[0], decomposed);
704 docolor=false;
705
706 }
707 if (GDL_DEBUG_PLSTREAM) fprintf(stderr,"draw_polyline()\n");
708 SizeT plotIndex=0;
709 bool line=false;
710 DLong local_psym=0;
711
712 if ( psym<0 )
713 {
714 line=true;
715 local_psym= -psym;
716 }
717 else if ( psym==0 )
718 {
719 line=true;
720 local_psym=psym;
721 }
722 else
723 {
724 local_psym=psym;
725 }
726
727 //usersym and other syms as well!
728 DFloat *userSymX, *userSymY;
729 DLong *userSymArrayDim;
730 //accelerate: define a localUsymArray where the computations that were previously in the loop are already done
731 DFloat *localUserSymX=NULL;
732 Guard<DFloat> guardlux;
733 DFloat *localUserSymY=NULL;
734 Guard<DFloat> guardluy;
735
736 bool useLocalPsymAccelerator=false;
737
738 //initialize symbol vertex list
739 static PLFLT xSym[49];
740 static PLFLT ySym[49];
741 DInt *do_fill;
742 bool *usersymhascolor;
743 bool *usersymhasthick;
744 DLong *usymColor;
745 DFloat *usymThick;
746 static DInt nofill=0;
747 if ( local_psym==8 )
748 {
749 GetUsym(&userSymArrayDim, &do_fill, &userSymX, &userSymY, &usersymhascolor, &usymColor , &usersymhasthick, &usymThick );
750 if ( *userSymArrayDim==0 )
751 {
752 ThrowGDLException("No user symbol defined.");
753 }
754 useLocalPsymAccelerator=true;
755 }
756 else if ( (local_psym>0&&local_psym<8))
757 {
758 do_fill=&nofill;
759 userSymArrayDim=&(syml[local_psym-1]);
760 useLocalPsymAccelerator=true;
761 switch(local_psym)
762 {
763 case 1:
764 userSymX=sym1x;
765 userSymY=sym1y;
766 break;
767 case 2:
768 userSymX=sym2x;
769 userSymY=sym2y;
770 break;
771 case 3:
772
773 //if device does not (bug?) draw single points (vector of 2 same coordinates as plplot does normally) then use special point psym=3.
774 if (GraphicsDevice::GetDevice()->DoesNotDrawSinglePoints()) {
775 userSymX = sym3xalt;
776 userSymY = sym3yalt;
777 userSymArrayDim=&syml_alt;
778 } else {
779 userSymX=sym3x;
780 userSymY=sym3y;
781 }
782 break;
783 case 4:
784 userSymX=sym4x;
785 userSymY=sym4y;
786 break;
787 case 5:
788 userSymX=sym5x;
789 userSymY=sym5y;
790 break;
791 case 6:
792 userSymX=sym6x;
793 userSymY=sym6y;
794 break;
795 case 7:
796 userSymX=sym7x;
797 userSymY=sym7y;
798 break;
799 }
800 }
801
802 if (useLocalPsymAccelerator) { //since userSymArrayDim is not defined
803 localUserSymX = (DFloat*) malloc(*userSymArrayDim * sizeof (DFloat));
804 guardlux.Reset(localUserSymX);
805 localUserSymY = (DFloat*) malloc(*userSymArrayDim * sizeof (DFloat));
806 guardluy.Reset(localUserSymY);
807 if (local_psym == 3 && GraphicsDevice::GetDevice()->DoesNotDrawSinglePoints()) {
808 for (int kk = 0; kk < *userSymArrayDim; kk++) {
809 localUserSymX[kk] = userSymX[kk] * a->getPsymConvX() / a->GetPlplotFudge() / a->getSymbolSize();
810 localUserSymY[kk] = userSymY[kk] * a->getPsymConvY() / a->GetPlplotFudge() / a->getSymbolSize();
811 }
812 } else {
813 for (int kk = 0; kk < *userSymArrayDim; kk++) {
814 localUserSymX[kk] = userSymX[kk] * a->getPsymConvX() / a->GetPlplotFudge();
815 localUserSymY[kk] = userSymY[kk] * a->getPsymConvY() / a->GetPlplotFudge();
816 }
817 }
818 }
819
820 DLong minEl=(xVal->N_Elements()<yVal->N_Elements())?
821 xVal->N_Elements():yVal->N_Elements();
822 // if scalar x
823 if ( xVal->N_Elements()==1&&xVal->Rank()==0 )
824 minEl=yVal->N_Elements();
825 // if scalar y
826 if ( yVal->N_Elements()==1&&yVal->Rank()==0 )
827 minEl=xVal->N_Elements();
828
829 // is one of the 2 "arrays" a singleton or not ?
830
831 PLFLT y, yMapBefore, y_ref;
832 int flag_y_const=0;
833 y_ref=static_cast<PLFLT>((*yVal)[0]);
834 if ( yVal->N_Elements()==1&&yVal->Rank()==0 ) flag_y_const=1;
835
836 PLFLT x, x1, xMapBefore, x_ref;
837 int flag_x_const=0;
838 x_ref=static_cast<PLFLT>((*xVal)[0]);
839 if ( xVal->N_Elements()==1&&xVal->Rank()==0 ) flag_x_const=1;
840
841 // AC 070601 we use a buffer to use the fast ->line method
842 // instead of the slow ->join one.
843 // 2 tricks:
844 // trick 1/ size of buffer is limited to 1e4 (compromize syze/speed) in order to be able to manage very
845 // large amount of data whitout duplicating all the arrays
846 // trick 2/ when we have a NaN or and Inf, we realize the plot, then reset.
847
848 int GDL_POLYLINE_BUFFSIZE=500000; // idl default seems to be more than 2e6 !!
849
850 if ( minEl<GDL_POLYLINE_BUFFSIZE ) GDL_POLYLINE_BUFFSIZE=append?minEl+1:minEl;
851 int i_buff=0;
852 PLFLT *x_buff=new PLFLT[GDL_POLYLINE_BUFFSIZE];
853 PLFLT *y_buff=new PLFLT[GDL_POLYLINE_BUFFSIZE];
854
855 bool isBad=FALSE;
856
857 for ( SizeT i=0; i<minEl; ++i ) {
858 isBad=FALSE;
859 if ( append ) //start with the old point
860 {
861 getLastPoint(a, x, y);
862 i--; //to get good counter afterwards
863 append=FALSE; //and stop appending after!
864 if ( xLog ) x=pow(10, x);
865 if ( yLog ) y=pow(10, y);
866 }
867 else
868 {
869 if ( !flag_x_const ) x=static_cast<PLFLT>((*xVal)[i]);
870 else x=x_ref;
871 if ( !flag_y_const ) y=static_cast<PLFLT>((*yVal)[i]);
872 else y=y_ref;
873 }
874
875 //note: here y is in minVal maxVal
876 if ( doMinMax ) isBad=((y<minVal)||(y>maxVal));
877 if ( xLog ) x=log10(x);
878 if ( yLog ) y=log10(y);
879 isBad=(isBad||!isfinite(x)|| !isfinite(y)||isnan(x)||isnan(y));
880 if ( isBad )
881 {
882 if ( i_buff>0 )
883 {
884 if ( line )
885 {
886 if (docolor) for (SizeT jj=0; jj< i_buff-1 ; ++jj)
887 {
888 a->Color ( ( *color )[plotIndex%color->N_Elements ( )], decomposed);
889 a->line(2, &(x_buff[jj]), &(y_buff[jj]));
890 plotIndex++;
891 }
892 else a->line(i_buff, x_buff, y_buff);
893 }
894 if (local_psym>0&&local_psym<8)
895 {
896 DLong oldStyl=gdlGetCurrentStyle();
897 a->styl(0, NULL, NULL); //symbols drawn in continuous lines
898 for ( int j=0; j<i_buff; ++j )
899 {
900 for ( int kk=0; kk < *userSymArrayDim; kk++ )
901 {
902 xSym[kk]=x_buff[j]+localUserSymX[kk];
903 ySym[kk]=y_buff[j]+localUserSymY[kk];
904 }
905 if (docolor)
906 {
907 a->Color ( ( *color )[plotIndex%color->N_Elements ( )], decomposed);
908 plotIndex++;
909 }
910 if ( *do_fill==1 )
911 {
912 a->fill(*userSymArrayDim, xSym, ySym);
913 }
914 else
915 {
916 a->line(*userSymArrayDim, xSym, ySym);
917 }
918 }
919 gdlLineStyle(a,oldStyl);
920 }
921 else if ( local_psym==8 )
922 {
923 DLong oldStyl=gdlGetCurrentStyle();
924 a->styl(0, NULL, NULL); //symbols drawn in continuous lines
925 if (*usersymhascolor)
926 {
927 a->Color(*usymColor, decomposed);
928 }
929 if (*usersymhasthick)
930 {
931 a->Thick(*usymThick);
932 }
933 for (int j = 0; j < i_buff; ++j)
934 {
935 for (int kk = 0; kk < *userSymArrayDim; kk++)
936 {
937 xSym[kk] = x_buff[j] + localUserSymX[kk];
938 ySym[kk] = y_buff[j] + localUserSymY[kk];
939 }
940 if (*do_fill == 1)
941 {
942 a->fill(*userSymArrayDim, xSym, ySym);
943 } else
944 {
945 a->line(*userSymArrayDim, xSym, ySym);
946 }
947 }
948 gdlLineStyle(a,oldStyl);
949 }
950 else if ( local_psym==10 )
951 {
952 ac_histo(a, i_buff, x_buff, y_buff, xLog);
953 }
954 i_buff=0;
955 }
956 continue;
957 }
958
959 x_buff[i_buff]=x;
960 y_buff[i_buff]=y;
961 i_buff=i_buff+1;
962
963 // cout << "nbuf: " << i << " " << i_buff << " "<< n_buff_max-1 << " " << minEl-1 << endl;
964
965 if ( (i_buff==GDL_POLYLINE_BUFFSIZE)||((i==minEl-1)&& !append)||((i==minEl)&&append) )
966 {
967 if ( line )
968 {
969 if (docolor) for (SizeT jj=0; jj< i_buff-1 ; ++jj)
970 {
971 a->Color ( ( *color )[plotIndex%color->N_Elements ( )], decomposed);
972 a->line(2, &(x_buff[jj]), &(y_buff[jj]));
973 plotIndex++;
974 }
975 else a->line(i_buff, x_buff, y_buff);
976 }
977 if ( local_psym>0&&local_psym<8 )
978 {
979 DLong oldStyl=gdlGetCurrentStyle();
980 a->styl(0, NULL, NULL); //symbols drawn in continuous lines
981 for ( int j=0; j<i_buff; ++j )
982 {
983 for ( int kk=0; kk < *userSymArrayDim; kk++ )
984 {
985 xSym[kk]=x_buff[j]+localUserSymX[kk];
986 ySym[kk]=y_buff[j]+localUserSymY[kk];
987 }
988 if (docolor)
989 {
990 a->Color ( ( *color )[plotIndex%color->N_Elements ( )], decomposed);
991 plotIndex++;
992 }
993 if ( *do_fill==1 )
994 {
995 a->fill(*userSymArrayDim, xSym, ySym);
996 }
997 else
998 {
999 a->line(*userSymArrayDim, xSym, ySym);
1000 }
1001 }
1002 gdlLineStyle(a,oldStyl);
1003 }
1004 else if ( local_psym==8 )
1005 {
1006 DLong oldStyl=gdlGetCurrentStyle();
1007 a->styl(0, NULL, NULL); //symbols drawn in continuous lines
1008 if (*usersymhascolor)
1009 {
1010 a->Color(*usymColor, decomposed);
1011 }
1012 if (*usersymhasthick)
1013 {
1014 a->Thick(*usymThick);
1015 }
1016 for (int j = 0; j < i_buff; ++j)
1017 {
1018 for (int kk = 0; kk < *userSymArrayDim; kk++)
1019 {
1020 xSym[kk] = x_buff[j] + localUserSymX[kk];
1021 ySym[kk] = y_buff[j] + localUserSymY[kk];
1022 }
1023 if (*do_fill == 1)
1024 {
1025 a->fill(*userSymArrayDim, xSym, ySym);
1026 } else
1027 {
1028 a->line(*userSymArrayDim, xSym, ySym);
1029 }
1030 }
1031 gdlLineStyle(a,oldStyl);
1032 }
1033 else if ( local_psym==10 )
1034 {
1035 ac_histo(a, i_buff, x_buff, y_buff, xLog);
1036 }
1037
1038 // we must recopy the last point since the line must continue (tested via small buffer ...)
1039 x_buff[0]=x_buff[i_buff-1];
1040 y_buff[0]=y_buff[i_buff-1];
1041 i_buff=1;
1042 }
1043 }
1044
1045 delete[] x_buff;
1046 delete[] y_buff;
1047 //save last point
1048 saveLastPoint(a, x, y);
1049 }
1050
1051
1052 //BACKGROUND COLOR
1053
1054
1055
1056 //Very special usage only in plotting surface
gdlSetGraphicsPenColorToBackground(GDLGStream * a)1057 void gdlSetGraphicsPenColorToBackground(GDLGStream *a)
1058 {
1059 a->plstream::col0( 0);
1060 }
1061
1062 //COLOR
1063
1064 // helper for NOERASE
1065
1066 //PSYM
1067
1068 //SYMSIZE
1069
1070 //CHARSIZE
1071
1072 //THICK
1073
1074 //crange to struct
1075
gdlStoreAxisCRANGE(int axisId,DDouble Start,DDouble End,bool log)1076 void gdlStoreAxisCRANGE(int axisId, DDouble Start, DDouble End, bool log)
1077 {
1078 DStructGDL* Struct=NULL;
1079 if ( axisId==XAXIS ) Struct=SysVar::X();
1080 if ( axisId==YAXIS ) Struct=SysVar::Y();
1081 if ( axisId==ZAXIS ) Struct=SysVar::Z();
1082 if ( Struct!=NULL )
1083 {
1084 int debug=0;
1085 if ( debug ) cout<<"Set :"<<Start<<" "<<End<<endl;
1086
1087 unsigned crangeTag=Struct->Desc()->TagIndex("CRANGE");
1088 if ( log )
1089 {
1090 (*static_cast<DDoubleGDL*>(Struct->GetTag(crangeTag, 0)))[0]=log10(Start);
1091 (*static_cast<DDoubleGDL*>(Struct->GetTag(crangeTag, 0)))[1]=log10(End);
1092 if ( debug ) cout<<"set log"<<Start<<" "<<End<<endl;
1093 }
1094 else
1095 {
1096 (*static_cast<DDoubleGDL*>(Struct->GetTag(crangeTag, 0)))[0]=Start;
1097 (*static_cast<DDoubleGDL*>(Struct->GetTag(crangeTag, 0)))[1]=End;
1098 }
1099 }
1100 }
1101
1102 //CRANGE from struct
1103
gdlGetCurrentAxisRange(int axisId,DDouble & Start,DDouble & End,bool checkMapset)1104 void gdlGetCurrentAxisRange(int axisId, DDouble &Start, DDouble &End, bool checkMapset)
1105 {
1106 DStructGDL* Struct=NULL;
1107 if ( axisId==XAXIS ) Struct=SysVar::X();
1108 if ( axisId==YAXIS ) Struct=SysVar::Y();
1109 if ( axisId==ZAXIS ) Struct=SysVar::Z();
1110 Start=0;
1111 End=0;
1112 if ( Struct!=NULL )
1113 {
1114 int debug=0;
1115 if ( debug ) cout<<"Get :"<<Start<<" "<<End<<endl;
1116 bool isProj;
1117 get_mapset(isProj);
1118 if (checkMapset && isProj && axisId!=ZAXIS) {
1119 DStructGDL* mapStruct=SysVar::Map(); //MUST NOT BE STATIC, due to .reset
1120 static unsigned uvboxTag=mapStruct->Desc()->TagIndex("UV_BOX");
1121 static DDoubleGDL *uvbox;
1122 uvbox=static_cast<DDoubleGDL*>(mapStruct->GetTag(uvboxTag, 0));
1123 if (axisId==XAXIS) {
1124 Start=(*uvbox)[0];
1125 End=(*uvbox)[2];
1126 } else {
1127 Start=(*uvbox)[1];
1128 End=(*uvbox)[3];
1129 }
1130 } else {
1131 static unsigned crangeTag=Struct->Desc()->TagIndex("CRANGE");
1132 Start=(*static_cast<DDoubleGDL*>(Struct->GetTag(crangeTag, 0)))[0];
1133 End=(*static_cast<DDoubleGDL*>(Struct->GetTag(crangeTag, 0)))[1];
1134
1135 static unsigned typeTag=Struct->Desc()->TagIndex("TYPE");
1136 if ( (*static_cast<DLongGDL*>(Struct->GetTag(typeTag, 0)))[0]==1 )
1137 {
1138 Start=pow(10., Start);
1139 End=pow(10., End);
1140 if ( debug ) cout<<"Get log :"<<Start<<" "<<End<<endl;
1141 }
1142 }
1143 }
1144 }
1145
gdlGetCurrentAxisWindow(int axisId,DDouble & wStart,DDouble & wEnd)1146 void gdlGetCurrentAxisWindow(int axisId, DDouble &wStart, DDouble &wEnd)
1147 {
1148 DStructGDL* Struct=NULL;
1149 if ( axisId==XAXIS ) Struct=SysVar::X();
1150 if ( axisId==YAXIS ) Struct=SysVar::Y();
1151 if ( axisId==ZAXIS ) Struct=SysVar::Z();
1152 wStart=0;
1153 wEnd=0;
1154 if ( Struct!=NULL )
1155 {
1156 static unsigned windowTag=Struct->Desc()->TagIndex("WINDOW");
1157 wStart=(*static_cast<DFloatGDL*>(Struct->GetTag(windowTag, 0)))[0];
1158 wEnd=(*static_cast<DFloatGDL*>(Struct->GetTag(windowTag, 0)))[1];
1159 }
1160 }
gdlStoreSC()1161 void gdlStoreSC() {
1162 //save corresponding SCxx values:
1163 DStructGDL* pStruct = SysVar::P(); //MUST NOT BE STATIC, due to .reset
1164 static unsigned positionTag = pStruct->Desc()->TagIndex("POSITION");
1165 DFloat* position = &(*static_cast<DFloatGDL*> (pStruct->GetTag(positionTag, 0)))[0];
1166 DStructGDL* dStruct = SysVar::D(); //MUST NOT BE STATIC, due to .reset
1167 static unsigned dxvsizeTag = dStruct->Desc()->TagIndex("X_VSIZE");
1168 static unsigned dyvsizeTag = dStruct->Desc()->TagIndex("Y_VSIZE");
1169
1170 DLong x_vsize = (*static_cast<DLongGDL*> (dStruct->GetTag(dxvsizeTag, 0)))[0];
1171 DLong y_vsize = (*static_cast<DLongGDL*> (dStruct->GetTag(dyvsizeTag, 0)))[0];
1172 DFloat* sc = SysVar::GetSC();
1173 DStructGDL* xStruct = SysVar::X();
1174 static unsigned xwindowTag = xStruct->Desc()->TagIndex("WINDOW");
1175 DStructGDL* yStruct = SysVar::Y();
1176 static unsigned ywindowTag = yStruct->Desc()->TagIndex("WINDOW");
1177 DFloat* xwindow=&(*static_cast<DFloatGDL*> (xStruct->GetTag(xwindowTag, 0)))[0];
1178 DFloat* ywindow=&(*static_cast<DFloatGDL*> (yStruct->GetTag(ywindowTag, 0)))[0];
1179 sc[0] = (position[2] != 0) ? position[0] * x_vsize : xwindow[0]*x_vsize;
1180 sc[1] = (position[2] != 0) ? position[2] * x_vsize : xwindow[1]*x_vsize;
1181 sc[2] = (position[2] != 0) ? position[1] * y_vsize : ywindow[0]*y_vsize;
1182 sc[3] = (position[2] != 0) ? position[3] * y_vsize : ywindow[1]*y_vsize;
1183
1184 }
1185 //Stores [XYZ].WINDOW, .REGION and .S
gdlStoreAxisSandWINDOW(GDLGStream * actStream,int axisId,DDouble Start,DDouble End,bool log)1186 void gdlStoreAxisSandWINDOW(GDLGStream* actStream, int axisId, DDouble Start, DDouble End, bool log)
1187 {
1188 PLFLT p_xmin, p_xmax, p_ymin, p_ymax, norm_min, norm_max, charDim;
1189 bool savesc=false;
1190 actStream->gvpd(p_xmin, p_xmax, p_ymin, p_ymax); //viewport normalized coords
1191 DStructGDL* Struct=NULL;
1192 if ( axisId==XAXIS ) {Struct=SysVar::X(); savesc=true; norm_min=p_xmin; norm_max=p_xmax; charDim=actStream->nCharLength();}
1193 if ( axisId==YAXIS ) {Struct=SysVar::Y(); savesc=true; norm_min=p_ymin; norm_max=p_ymax; charDim=actStream->nCharHeight();}
1194 if ( axisId==ZAXIS ) {Struct=SysVar::Z(); norm_min=0; norm_max=1; charDim=actStream->nCharLength();}
1195 if ( Struct!=NULL )
1196 {
1197 unsigned marginTag=Struct->Desc()->TagIndex("MARGIN");
1198 DFloat m1=(*static_cast<DFloatGDL*>(Struct->GetTag(marginTag, 0)))[0];
1199 DFloat m2=(*static_cast<DFloatGDL*>(Struct->GetTag(marginTag, 0)))[1];
1200 static unsigned regionTag=Struct->Desc()->TagIndex("REGION");
1201 (*static_cast<DFloatGDL*>(Struct->GetTag(regionTag, 0)))[0]=max(0.0,norm_min-m1*charDim);
1202 (*static_cast<DFloatGDL*>(Struct->GetTag(regionTag, 0)))[1]=min(1.0,norm_max+m2*charDim);
1203
1204 // if ( log ) {Start=log10(Start); End=log10(End);}
1205 static unsigned windowTag=Struct->Desc()->TagIndex("WINDOW");
1206 (*static_cast<DFloatGDL*>(Struct->GetTag(windowTag, 0)))[0]=norm_min;
1207 (*static_cast<DFloatGDL*>(Struct->GetTag(windowTag, 0)))[1]=norm_max;
1208
1209 static unsigned sTag=Struct->Desc()->TagIndex("S");
1210 (*static_cast<DDoubleGDL*>(Struct->GetTag(sTag, 0)))[0]=
1211 (norm_min*End-norm_max*Start)/(End-Start);
1212 (*static_cast<DDoubleGDL*>(Struct->GetTag(sTag, 0)))[1]=
1213 (norm_max-norm_min)/(End-Start);
1214 if (savesc) gdlStoreSC();
1215 }
1216 }
1217
gdlStoreCLIP(DLongGDL * clipBox)1218 void gdlStoreCLIP(DLongGDL* clipBox)
1219 {
1220 DStructGDL* pStruct=SysVar::P(); //MUST NOT BE STATIC, due to .reset
1221 int i;
1222 static unsigned clipTag=pStruct->Desc()->TagIndex("CLIP");
1223 for ( i=0; i<clipBox->N_Elements(); ++i ) (*static_cast<DLongGDL*>(pStruct->GetTag(clipTag, 0)))[i]=(*clipBox)[i];
1224 }
1225
gdlGetAxisType(int axisId,bool & log)1226 void gdlGetAxisType(int axisId, bool &log)
1227 {
1228 DStructGDL* Struct;
1229 if ( axisId==XAXIS ) Struct=SysVar::X();
1230 if ( axisId==YAXIS ) Struct=SysVar::Y();
1231 if ( axisId==ZAXIS ) Struct=SysVar::Z();
1232 if ( Struct!=NULL )
1233 {
1234 static unsigned typeTag=Struct->Desc()->TagIndex("TYPE");
1235 if ( (*static_cast<DLongGDL*>(Struct->GetTag(typeTag, 0)))[0]==1 )
1236 log=true;
1237 else
1238 log=false;
1239 }
1240 }
1241
get_mapset(bool & mapset)1242 void get_mapset(bool &mapset)
1243 {
1244 DStructGDL* Struct=SysVar::X(); //MUST NOT BE STATIC, due to .reset
1245 if ( Struct!=NULL )
1246 {
1247 static unsigned typeTag=Struct->Desc()->TagIndex("TYPE");
1248
1249 if ( (*static_cast<DLongGDL*>(Struct->GetTag(typeTag, 0)))[0]==3 )
1250 mapset=true;
1251 else
1252 mapset=false;
1253 }
1254 }
1255
set_mapset(bool mapset)1256 void set_mapset(bool mapset)
1257 {
1258 DStructGDL* Struct=SysVar::X(); //MUST NOT BE STATIC, due to .reset
1259 if ( Struct!=NULL )
1260 {
1261 static unsigned typeTag=Struct->Desc()->TagIndex("TYPE");
1262 (*static_cast<DLongGDL*>(Struct->GetTag(typeTag, 0)))[0]=(mapset)?3:0;
1263 }
1264 }
1265
1266
1267 //axis type (log..)
1268
gdlStoreAxisType(int axisId,bool Type)1269 void gdlStoreAxisType(int axisId, bool Type)
1270 {
1271 DStructGDL* Struct=NULL;
1272 if ( axisId==XAXIS ) Struct=SysVar::X();
1273 if ( axisId==YAXIS ) Struct=SysVar::Y();
1274 if ( axisId==ZAXIS ) Struct=SysVar::Z();
1275 if ( Struct!=NULL )
1276 {
1277 static unsigned typeTag=Struct->Desc()->TagIndex("TYPE");
1278 (*static_cast<DLongGDL*>(Struct->GetTag(typeTag, 0)))[0]=Type;
1279 }
1280 }
1281
doOurOwnFormat(PLINT axisNotUsed,PLFLT value,char * label,PLINT length,PLPointer data)1282 void doOurOwnFormat(PLINT axisNotUsed, PLFLT value, char *label, PLINT length, PLPointer data)
1283 {
1284 struct GDL_TICKDATA *ptr = (GDL_TICKDATA* )data;
1285 //was:
1286 // static string normalfmt[7]={"%1.0fx10#u%d#d","%2.1fx10#u%d#d","%3.2fx10#u%d#d","%4.2fx10#u%d#d","%5.4fx10#u%d#d","%6.5fx10#u%d#d","%7.6fx10#u%d#d"};
1287 // static string specialfmt="10#u%d#d";
1288 // static string specialfmtlog="10#u%s#d";
1289
1290 //we need !3x!X to insure the x sign is always written in single roman.
1291 static string normalfmt[7]={"%1.0f!3x!X10!E%d!N","%2.1f!3x!X10!E%d!N","%3.2f!3x!X10!E%d!N","%4.2f!3x!X10!E%d!N","%5.4f!3x!X10!E%d!N","%6.5f!3x!X10!E%d!N","%7.6f!3x!X10!E%d!N"};
1292 static string specialfmt="10!E%d!N";
1293 static string specialfmtlog="10!E%s!N";
1294 PLFLT z;
1295 int ns;
1296 char *i;
1297 int sgn=(value<0)?-1:1;
1298 //special cases, since plplot gives approximate zero values, not strict zeros.
1299 if (!(ptr->isLog) && (sgn*value<ptr->axisrange*1e-6))
1300 {
1301 snprintf(label, length, "0");
1302 return;
1303 }
1304 //in log, plplot gives correctly rounded "integer" values but 10^0 needs a bit of help.
1305 if ((ptr->isLog) && (sgn*value<1e-6)) //i.e. 0
1306 {
1307 snprintf(label, length, "1");
1308 return;
1309 }
1310
1311 int e=floor(log10(value*sgn));
1312 char *test=(char*)calloc(2*length, sizeof(char)); //be safe
1313 if (!isfinite(log10(value*sgn))||(e<4 && e>-4))
1314 {
1315 snprintf(test, length, "%f",value);
1316 ns=strlen(test);
1317 i=strrchr (test,'0');
1318 while (i==(test+ns-1)) //remove trailing zeros...
1319 {
1320 *i='\0';
1321 i=strrchr(test,'0');
1322 ns--;
1323 }
1324 i=strrchr(test,'.'); //remove trailing '.'
1325 if (i==(test+ns-1)) {*i='\0'; ns--;}
1326 if (ptr->isLog) {
1327 snprintf( label, length, specialfmtlog.c_str(),test);
1328 }
1329 else
1330 {
1331 strcpy(label, test);
1332 }
1333 }
1334 else
1335 {
1336 z=value*sgn/pow(10.,e);
1337 snprintf(test,20,"%7.6f",z);
1338 ns=strlen(test);
1339 i=strrchr(test,'0');
1340 while (i==(test+ns-1))
1341 {
1342 *i='\0';
1343 i=strrchr(test,'0');
1344 ns--;
1345 }
1346 ns-=2;ns=(ns>6)?6:ns;
1347 if (floor(sgn*z)==1 && ns==0) {
1348 snprintf( label, length, specialfmt.c_str(),e);
1349 } else {
1350 snprintf( label, length, normalfmt[ns].c_str(),sgn*z,e);
1351 }
1352 }
1353 free(test);
1354 }
doFormatAxisValue(DDouble value,string & label)1355 void doFormatAxisValue(DDouble value, string &label)
1356 {
1357 static string normalfmt[7]={"%1.0fx10^%d","%2.1fx10^%d","%3.2fx10^%d","%4.2fx10^%d","%5.4fx10^%d","%6.5fx10^%d","%7.6fx10^%d"};
1358 static string specialfmt="10^%d";
1359 static const int length=20;
1360 PLFLT z;
1361 int ns;
1362 char *i;
1363 int sgn=(value<0)?-1:1;
1364 //special cases, since plplot gives approximate zero values, not strict zeros.
1365 if (sgn*value< std::numeric_limits<DDouble>::min())
1366 {
1367 label="0";
1368 return;
1369 }
1370
1371 int e=floor(log10(value*sgn));
1372 char *test=(char*)calloc(2*length, sizeof(char)); //be safe
1373 if (!isfinite(e)||(e<4 && e>-4))
1374 {
1375 snprintf(test, length, "%f",value);
1376 ns=strlen(test);
1377 i=strrchr (test,'0');
1378 while (i==(test+ns-1)) //remove trailing zeros...
1379 {
1380 *i='\0';
1381 i=strrchr(test,'0');
1382 ns--;
1383 }
1384 i=strrchr(test,'.'); //remove trailing '.'
1385 if (i==(test+ns-1)) {*i='\0'; ns--;}
1386 }
1387 else
1388 {
1389 z=value*sgn/pow(10.,e);
1390 snprintf(test,20,"%7.6f",z);
1391 ns=strlen(test);
1392 i=strrchr(test,'0');
1393 while (i==(test+ns-1))
1394 {
1395 *i='\0';
1396 i=strrchr(test,'0');
1397 ns--;
1398 }
1399 ns-=2;ns=(ns>6)?6:ns;
1400 if (floor(sgn*z)==1 && ns==0)
1401 snprintf( test, length, specialfmt.c_str(),e);
1402 else
1403 snprintf( test, length, normalfmt[ns].c_str(),sgn*z,e);
1404 }
1405 label=test;
1406 free(test);
1407 }
1408
format_axis_values(EnvT * e)1409 BaseGDL* format_axis_values(EnvT *e){
1410 DDoubleGDL* p0D = e->GetParAs<DDoubleGDL>( 0);
1411 DStringGDL* res = new DStringGDL( p0D->Dim(), BaseGDL::NOZERO);
1412 SizeT nEl = p0D->N_Elements();
1413 for( SizeT i=0; i<nEl; ++i)
1414 {
1415 doFormatAxisValue((*p0D)[i], (*res)[i]);
1416 }
1417 return res;
1418 }
1419 //just a wrapper for doOurOwnFormat() adding general font code translation.
gdlSimpleAxisTickFunc(PLINT axis,PLFLT value,char * label,PLINT length,PLPointer data)1420 void gdlSimpleAxisTickFunc(PLINT axis, PLFLT value, char *label, PLINT length, PLPointer data)
1421 {
1422 struct GDL_TICKDATA *ptr = (GDL_TICKDATA* )data;
1423 doOurOwnFormat(axis, value, label, length, data);
1424 //translate format codes (as in mtex).
1425 double nchars;
1426 std::string out = ptr->a->TranslateFormatCodes(label, &nchars);
1427 ptr->nchars=max(ptr->nchars,nchars);
1428 strcpy(label,out.c_str());
1429 }
1430
gdlMultiAxisTickFunc(PLINT axis,PLFLT value,char * label,PLINT length,PLPointer multiaxisdata)1431 void gdlMultiAxisTickFunc(PLINT axis, PLFLT value, char *label, PLINT length, PLPointer multiaxisdata)
1432 {
1433 PLINT axisNotUsed=0; //to indicate that effectively the axis number is not (yet?) used in some calls
1434 static GDL_TICKDATA tdata;
1435 static SizeT internalIndex=0;
1436 static DLong lastMultiAxisLevel=0;
1437 static string theMonth[12]={"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"};
1438 PLINT Month, Day , Year , Hour , Minute, dow, cap;
1439 PLFLT Second;
1440 struct GDL_MULTIAXISTICKDATA *ptr = (GDL_MULTIAXISTICKDATA* )multiaxisdata;
1441 tdata.a=ptr->a;
1442 tdata.isLog=ptr->isLog;
1443 if (ptr->reset) {
1444 internalIndex=0; //reset index each time a new axis command is issued.
1445 ptr->reset=false;
1446 }
1447 if (ptr->counter != lastMultiAxisLevel)
1448 {
1449 lastMultiAxisLevel=ptr->counter;
1450 internalIndex=0; //reset index each time sub-axis changes
1451 }
1452
1453 if (ptr->what==GDL_TICKFORMAT || (ptr->what==GDL_TICKFORMAT_AND_UNITS && ptr->counter < ptr->nTickFormat) )
1454 {
1455 if (ptr->counter > ptr->nTickFormat-1)
1456 {
1457 doOurOwnFormat(axisNotUsed, value, label, length, &tdata);
1458 }
1459 else
1460 { //must pass the value, not the log, to the formatter?
1461 DDouble v=value;
1462 if (tdata.isLog) v=pow(10.,v);
1463 if (((*ptr->TickFormat)[ptr->counter]).substr(0,1) == "(")
1464 { //internal format, call internal func "STRING"
1465 EnvT *e=ptr->e;
1466 static int stringIx = LibFunIx("STRING");
1467 assert( stringIx >= 0);
1468 EnvT* newEnv= new EnvT(e, libFunList[stringIx], NULL);
1469 Guard<EnvT> guard( newEnv);
1470 // add parameters
1471 newEnv->SetNextPar( new DDoubleGDL(v));
1472 newEnv->SetNextPar( new DStringGDL(((*ptr->TickFormat)[ptr->counter]).c_str()));
1473 // make the call
1474 BaseGDL* res = static_cast<DLibFun*>(newEnv->GetPro())->Fun()(newEnv);
1475 strncpy(label,(*static_cast<DStringGDL*>(res))[0].c_str(),1000);
1476 }
1477 else // external function: if tickunits not specified, pass Axis (int), Index(int),Value(Double)
1478 // else pass also Level(int)
1479 // Thanks to Marc for code snippet!
1480 // NOTE: this encompasses the 'LABEL_DATE' format, an existing procedure in the IDL library.
1481 {
1482 EnvT *e=ptr->e;
1483 DString callF=(*ptr->TickFormat)[ptr->counter];
1484 // this is a function name -> convert to UPPERCASE
1485 callF = StrUpCase( callF);
1486 // Search in user proc and function
1487 SizeT funIx = GDLInterpreter::GetFunIx( callF);
1488
1489 EnvUDT* newEnv = new EnvUDT( e->CallingNode(), funList[ funIx], (DObjGDL**)NULL);
1490 Guard< EnvUDT> guard( newEnv);
1491 // add parameters
1492 newEnv->SetNextPar( new DLongGDL(axis-1)); //axis in PLPLOT starts at 1, it starts at 0 in IDL
1493 newEnv->SetNextPar( new DLongGDL(internalIndex)); //index
1494 newEnv->SetNextPar( new DDoubleGDL(v)); //value
1495 if (ptr->what==GDL_TICKFORMAT_AND_UNITS) newEnv->SetNextPar( new DLongGDL(ptr->counter)); //level
1496 // guard *before* pushing new env
1497 StackGuard<EnvStackT> guard1 ( e->Interpreter()->CallStack());
1498 e->Interpreter()->CallStack().push_back(newEnv);
1499 guard.release();
1500
1501 BaseGDL* retValGDL = e->Interpreter()->call_fun(static_cast<DSubUD*>(newEnv->GetPro())->GetTree());
1502 // we are the owner of the returned value
1503 Guard<BaseGDL> retGuard( retValGDL);
1504 strncpy(label,(*static_cast<DStringGDL*>(retValGDL))[0].c_str(),1000);
1505 }
1506 }
1507 }
1508 else if (ptr->what==GDL_TICKUNITS || (ptr->what==GDL_TICKFORMAT_AND_UNITS && ptr->counter >= ptr->nTickFormat))
1509 {
1510 if (ptr->counter > ptr->nTickUnits-1 )
1511 {
1512 doOurOwnFormat(axisNotUsed, value, label, length, &tdata);
1513 }
1514 else
1515 {
1516 DString what=StrUpCase((*ptr->TickUnits)[ptr->counter]);
1517 if (what.length()<1) {
1518 doOurOwnFormat(axisNotUsed, value, label, length, &tdata);
1519 }
1520 else if (what.substr(0,7)=="NUMERIC") {
1521 doOurOwnFormat(axisNotUsed, value, label, length, &tdata);
1522 } else {
1523 j2ymdhms(value, Month , Day , Year , Hour , Minute, Second, dow, cap);
1524 int convcode=0;
1525 if (what.length()<1) convcode=7;
1526 else if (what.substr(0,4)=="YEAR") convcode=1;
1527 else if (what.substr(0,5)=="MONTH") convcode=2;
1528 else if (what.substr(0,3)=="DAY") convcode=3;
1529 else if (what.substr(0,4)=="HOUR") convcode=4;
1530 else if (what.substr(0,6)=="MINUTE") convcode=5;
1531 else if (what.substr(0,6)=="SECOND") convcode=6;
1532 else if (what.substr(0,4)=="TIME")
1533 {
1534 if(ptr->axisrange>=366) convcode=1;
1535 else if(ptr->axisrange>=32) convcode=2;
1536 else if(ptr->axisrange>=1.1) convcode=3;
1537 else if(ptr->axisrange*24>=1.1) convcode=4;
1538 else if(ptr->axisrange*24*60>=1.1) convcode=5;
1539 else convcode=6;
1540 } else convcode=7;
1541 switch(convcode){
1542 case 1:
1543 snprintf( label, length, "%d", Year);
1544 break;
1545 case 2:
1546 snprintf( label, length, "%s", theMonth[Month].c_str());
1547 break;
1548 case 3:
1549 snprintf( label, length, "%d", Day);
1550 break;
1551 case 4:
1552 snprintf( label, length, "%02d", Hour);
1553 break;
1554 case 5:
1555 snprintf( label, length, "%02d", Minute);
1556 break;
1557 case 6:
1558 snprintf( label, length, "%05.2f", Second);
1559 break;
1560 case 7:
1561 doOurOwnFormat(axisNotUsed, value, label, length, &tdata);
1562 break;
1563 }
1564
1565 }
1566 }
1567 }
1568 //translate format codes (as in mtex).
1569 double nchars;
1570 std::string out = ptr->a->TranslateFormatCodes(label, &nchars);
1571 ptr->nchars=max(ptr->nchars,nchars);
1572 strcpy(label,out.c_str());
1573 internalIndex++;
1574 }
1575
gdlSingleAxisTickNamedFunc(PLINT axisNotUsed,PLFLT value,char * label,PLINT length,PLPointer data)1576 void gdlSingleAxisTickNamedFunc( PLINT axisNotUsed, PLFLT value, char *label, PLINT length, PLPointer data)
1577 {
1578 static GDL_TICKDATA tdata;
1579 struct GDL_TICKNAMEDATA *ptr = (GDL_TICKNAMEDATA* )data;
1580 tdata.isLog=ptr->isLog;
1581 tdata.axisrange=ptr->axisrange;
1582 if (ptr->counter > ptr->nTickName-1)
1583 {
1584 doOurOwnFormat(axisNotUsed, value, label, length, &tdata);
1585 }
1586 else
1587 {
1588 snprintf( label, length, "%s", ((*ptr->TickName)[ptr->counter]).c_str() );
1589 }
1590 //translate format codes (as in mtex).
1591 double nchars;
1592 std::string out = ptr->a->TranslateFormatCodes(label, &nchars);
1593 ptr->nchars=max(ptr->nchars,nchars);
1594 strcpy(label,out.c_str());
1595 ptr->counter++;
1596 }
1597
T3Denabled()1598 bool T3Denabled()
1599 {
1600 DStructGDL* pStruct=SysVar::P(); //MUST NOT BE STATIC, due to .reset
1601 DLong ok4t3d=(*static_cast<DLongGDL*>(pStruct->GetTag(pStruct->Desc()->TagIndex("T3D"), 0)))[0];
1602 if (ok4t3d==0) return false; else return true;
1603 }
1604
usersym(EnvT * e)1605 void usersym(EnvT *e)
1606 {
1607 DFloatGDL *xyVal, *xVal, *yVal;
1608 Guard<BaseGDL> p0_guard;
1609 DLong n;
1610 DInt do_fill;
1611 bool do_color;
1612 bool do_thick;
1613 DFloat thethick;
1614 DLong thecolor;
1615 DFloat *x, *y;
1616 SizeT nParam=e->NParam();
1617 if (nParam==0) e->Throw("Incorrect number of arguments.");
1618 if ( nParam==1 )
1619 {
1620 BaseGDL* p0=e->GetNumericArrayParDefined(0)->Transpose(NULL); //hence [49,2]
1621
1622 xyVal=static_cast<DFloatGDL*>
1623 (p0->Convert2(GDL_FLOAT, BaseGDL::COPY));
1624 p0_guard.Reset(p0); // delete upon exit
1625
1626 if ( xyVal->Rank()!=2||xyVal->Dim(1)!=2 )
1627 e->Throw(e->GetParString(0)+" must be a 2-dim array of type [2,N] in this context.");
1628
1629 if ( xyVal->Dim(0)>49 )
1630 {
1631 e->Throw("Max array size for USERSYM is 49");
1632 }
1633 n=xyVal->Dim(0);
1634 // array is in the good order for direct C assignement
1635 x=&(*xyVal)[0];
1636 y=&(*xyVal)[n];
1637 }
1638 else
1639 {
1640 xVal=e->GetParAs< DFloatGDL>(0);
1641 if ( xVal->Rank()!=1 )
1642 e->Throw(e->GetParString(0)+" must be a 1D array in this context: ");
1643
1644 yVal=e->GetParAs< DFloatGDL>(1);
1645 if ( yVal->Rank()!=1 )
1646 e->Throw("Expression must be a 1D array in this context: "+e->GetParString(1));
1647
1648 if ( xVal->Dim(0)!=yVal->Dim(0) )
1649 {
1650 e->Throw("Arrays must have same size ");
1651 }
1652
1653 if ( xVal->Dim(0)>49 )
1654 {
1655 e->Throw("Max array size for USERSYM is 49");
1656 }
1657 n=xVal->Dim(0);
1658 x=&(*xVal)[0];
1659 y=&(*yVal)[0];
1660 }
1661 do_fill=0;
1662 static int FILLIx = e->KeywordIx("FILL");
1663 if ( e->KeywordSet(FILLIx) )
1664 {
1665 do_fill=1;
1666 }
1667 //IDL does not complain if color is undefined.
1668 do_color=false;
1669 thecolor=0;
1670 static int COLORIx = e->KeywordIx("COLOR");
1671 if ( e->KeywordPresent(COLORIx))
1672 {
1673 if (e->IfDefGetKWAs<DLongGDL>( COLORIx )) {
1674 e->AssureLongScalarKW(COLORIx, thecolor);
1675 do_color=true;
1676 }
1677 }
1678 //IDL does not complain if thick is undefined.
1679 do_thick=false;
1680 thethick=0;
1681 static int THICKIx = e->KeywordIx("THICK");
1682 if ( e->KeywordPresent(THICKIx))
1683 {
1684 if (e->IfDefGetKWAs<DFloatGDL>( THICKIx )) {
1685 e->AssureFloatScalarKW(THICKIx, thethick);
1686 do_thick=true;
1687 }
1688 }
1689 SetUsym(n, do_fill, x, y, do_color, thecolor, do_thick, thethick);
1690 }
1691 #ifdef USE_LIBPROJ
GDLgrProjectedPolygonPlot(GDLGStream * a,PROJTYPE ref,DStructGDL * map,DDoubleGDL * lons_donottouch,DDoubleGDL * lats_donottouch,bool isRadians,bool const doFill,DLongGDL * conn)1692 void GDLgrProjectedPolygonPlot( GDLGStream * a, PROJTYPE ref, DStructGDL* map,
1693 DDoubleGDL *lons_donottouch, DDoubleGDL *lats_donottouch, bool isRadians, bool const doFill, DLongGDL *conn ) {
1694 DDoubleGDL *lons,*lats;
1695 lons=lons_donottouch->Dup(); Guard<BaseGDL> lonsGuard( lons);
1696 lats=lats_donottouch->Dup(); Guard<BaseGDL> latsGuard( lats);
1697
1698 DStructGDL* localMap = map;
1699 if (localMap==NULL) localMap=SysVar::Map( );
1700 bool mapSet;
1701 get_mapset(mapSet); //if mapSet, output will be converted to normalized coordinates as this seems to be the way to do it.
1702 bool doConn = (conn != NULL);
1703 DLongGDL *gons, *lines;
1704 if (!isRadians) {
1705 SizeT nin = lons->N_Elements( );
1706 #pragma omp parallel if (nin >= CpuTPOOL_MIN_ELTS && (CpuTPOOL_MAX_ELTS == 0 || CpuTPOOL_MAX_ELTS <= nin))
1707 {
1708 #pragma omp for
1709 for ( OMPInt in = 0; in < nin; in++ ) { //pass in radians for gdlProjForward
1710 (*lons)[in] *= DEG_TO_RAD;
1711 (*lats)[in] *= DEG_TO_RAD;
1712 }
1713 }
1714 }
1715 DDoubleGDL *res = gdlProjForward( ref, localMap, lons, lats, conn, doConn, gons, doFill, lines, !doFill, false );
1716 SizeT nout = res->N_Elements( ) / 2;
1717 if (nout < 1) {GDLDelete(res); return;} //projection clipped totally these values.
1718 DDoubleGDL *res2 = static_cast<DDoubleGDL*> (static_cast<BaseGDL*> (res)->Transpose( NULL )); GDLDelete(res);
1719 int minpoly;
1720 if ( doFill ) {
1721 conn = gons;
1722 minpoly = 3;
1723 } else {
1724 conn = lines;
1725 minpoly = 2;
1726 }
1727 SizeT index = 0;
1728 SizeT size;
1729 SizeT start;
1730 while ( index < conn->N_Elements( ) ) {
1731 size = (*conn)[index];
1732 if ( size == 0 ) break; //cannot be negative!
1733 start = (*conn)[index + 1];
1734 if ( size >= minpoly ) {
1735 if ( doFill ) {
1736 a->fill( size, (PLFLT*) &((*res2)[start]), (PLFLT*) &((*res2)[start + nout]) );
1737 } else {
1738 a->line( size, (PLFLT*) &((*res2)[start]), (PLFLT*) &((*res2)[start + nout]) );
1739 }
1740 }
1741 index += (size + 1);
1742 }
1743 GDLDelete( res2 );
1744 if ( doFill ) GDLDelete( gons );
1745 else GDLDelete( lines );
1746 }
1747 #endif
1748 } // namespace
1749