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