home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 449a.lha / ListPlot / Csrc / plotting.c < prev   
Encoding:
C/C++ Source or Header  |  1990-12-08  |  24.9 KB  |  1,165 lines
  1. /*
  2.  * plotting.c
  3.  *
  4.  * plotting routines.
  5.  *
  6.  */
  7. #include <stdio.h>
  8. #include <errno.h>
  9. extern int errno;
  10. #include <assert.h>
  11. #include <string.h>
  12. #include <ctype.h>
  13. #include <signal.h>
  14. #include "datatypes.h"
  15.  
  16. void    ErrorExit();
  17.  
  18. #ifndef    min
  19. #    define    min(A,B)    A<B? A : B
  20. #endif
  21. #ifndef    max
  22. #    define    max(A,B)    A>B? A : B
  23. #endif
  24.  
  25. extern bool    GraphicsInProgress;
  26.  
  27. /* GoldenRatio = (1+Sqrt(5))/2    */
  28. #define    GOLDENRATIO    1.6180339887499
  29. #define ASPECTRATIO    1.0/GOLDENRATIO
  30. extern VALUE    V_AngularUnit;
  31. extern VALUE    V_AnnotationScale;
  32. extern VALUE    V_AspectRatio;
  33. extern VALUE    V_Boxed;
  34. extern VALUE    V_Domain;
  35. extern VALUE    V_Gridding;
  36. extern VALUE    V_LabelScale;
  37. extern VALUE    V_LineColor;
  38. extern VALUE    V_LineStyle;
  39. extern VALUE    V_Orientation;
  40. extern VALUE    V_Origin;
  41. extern VALUE    V_PlotColor;
  42. extern VALUE    V_PlotDevice;
  43. extern VALUE    V_PlotJoined;
  44. extern VALUE    V_PlotPoints;
  45. extern VALUE    V_PlotTitle;
  46. extern VALUE    V_PlotType;
  47. extern VALUE    V_PointScale;
  48. extern VALUE    V_PointSymbol;
  49. extern VALUE    V_PolarVariable;
  50. extern VALUE    V_Range;
  51. extern VALUE    V_SubPages;
  52. extern VALUE    V_SupplyAbscissa;
  53. extern VALUE    V_TitleScale;
  54. extern VALUE    V_UseInputFile;
  55. extern VALUE    V_Verbose;
  56. extern VALUE    V_ViewPort;
  57. extern VALUE    V_XLabel;
  58. extern VALUE    V_XTick;
  59. extern VALUE    V_YLabel;
  60. extern VALUE    V_YTick;
  61.  
  62. char    *DeviceTypes[] = {
  63.     "amiga",
  64.     "printer",
  65.     "iff",
  66.     "hp",
  67.     "aegis",
  68.     "postscript",
  69.     (char *)NULL
  70. };
  71.  
  72. typedef struct line_style {
  73.         int    ls_type;
  74.         int    ls_nelms;
  75.         int    *ls_space;
  76.         int    *ls_mark;
  77.     } LINESTYLE;
  78. LIST    LineStyleList;
  79. LINESTYLE    *LineStyles; int NLineStyles;
  80. int        *PointCodes, NPointCodes;
  81. int    space0 = 0, mark0 = 0, space1 = 1500, mark1 = 1500;
  82.  
  83. RECT    UserRect;
  84. RECT    ViewPort = { 0.1,0.1, 0.9,0.9 };
  85. RECT    WorldRect = {0.0,0.0, 1.0,1.0};
  86.  
  87.  
  88. void
  89. ListPlot(Data, NTuples, TupleSize)
  90. FLOAT    **Data;
  91. int    NTuples, TupleSize;
  92. {
  93.     void InitializeDevice();
  94.     void InitializePlottingEnv();
  95. #ifdef    ANSI_C
  96.     void Plot2D(FLOAT **, int , int );
  97.     void PolarPlot(FLOAT **, int , int );
  98. #else
  99.     void Plot2D();
  100.     void PolarPlot();
  101. #endif
  102.  
  103.     InitializeDevice();
  104.     GraphicsInProgress = TRUE;
  105.  
  106.     InitializePlottingEnv(Data, NTuples, TupleSize);
  107.  
  108.     if (Vstrcmp(V_PlotType,"polar")) {
  109.         PolarPlot(Data, NTuples, TupleSize);
  110.     } else {
  111.         Plot2D(Data, NTuples, TupleSize);
  112.     }
  113.     plend();
  114. }
  115.  
  116.  
  117. void
  118. InitializeDevice()
  119. {
  120.     register int    i;
  121.     register char    *DeviceStr;
  122.  
  123.  
  124.     plorient(Vstrcmp(V_Orientation, "portrait")? 1 : 0);
  125.     plselect(stdout);    /* write to stdout */
  126.         
  127.  
  128.     DeviceStr = V_PlotDevice.val_u.udt_string;
  129.     for (i=0; DeviceTypes[i]; i++) 
  130.         if (strcmp(DeviceStr, DeviceTypes[i]) == 0) {
  131.             plbeg(
  132.                 i+1,    /* device code    */
  133.                 (int)V_SubPages.val_u.udt_interval.int_lo, /* nx */
  134.                 (int)V_SubPages.val_u.udt_interval.int_hi /* ny */
  135.             );
  136.             return;
  137.         }
  138.     fprintf(stderr,"(InitializeDevice) Unknown device type: %s\n", DeviceStr);
  139.     ErrorExit();
  140. }
  141.  
  142. void
  143. InitializePlottingEnv(Data, NTuples, TupleSize)
  144. FLOAT    **Data;
  145. int    NTuples;
  146. int    TupleSize;
  147. {
  148.     int    i, j;
  149.     FLOAT    *Ptr;
  150.     void    InitLineStyles();
  151.     void    InitPointCodes();
  152.  
  153.     /* viewport    */
  154.     ViewPort = VtoRect(V_ViewPort);
  155.  
  156.     /* user window        */
  157.     UserRect.XMin = MAX_FLOAT;
  158.     UserRect.YMin = MAX_FLOAT;
  159.     UserRect.XMax = -MAX_FLOAT;
  160.     UserRect.YMax = -MAX_FLOAT;
  161.     
  162.     if (VtoBoolean(V_SupplyAbscissa)) {
  163.         UserRect.XMin = 0.0;
  164.         UserRect.XMax = (FLOAT)(NTuples-1);
  165.     } else {
  166.         for (i=0; i<NTuples; i++) {
  167.             if (Data[i][0] < UserRect.XMin)
  168.                 UserRect.XMin = Data[i][0];
  169.             if (Data[i][0] > UserRect.XMax)
  170.                 UserRect.XMax = Data[i][0];
  171.         }
  172.     }
  173.     for (i=0; i<NTuples; i++) 
  174.         for (j=TupleSize-1,Ptr= &Data[i][1]; j>0; --j,Ptr++) {
  175.             if (*Ptr < UserRect.YMin)
  176.                 UserRect.YMin = *Ptr;
  177.             if (*Ptr > UserRect.YMax)
  178.                 UserRect.YMax = *Ptr;
  179.         }
  180.         
  181.     
  182.     /* font scaling        */
  183.  
  184.     /* plot type        */
  185.  
  186.     /* grid            */
  187.  
  188.     /* line styles        */
  189.     InitLineStyles();
  190.     
  191.  
  192.     /* line colors        */
  193.  
  194.     /* point symbols    */
  195.     InitPointCodes();
  196.  
  197.     /* window type        */
  198.  
  199.     /* X label        */
  200.  
  201.     /* Y label        */
  202.  
  203.     /* Plot title        */
  204.  
  205.     /* aspect ratio        */
  206. }
  207.  
  208. #define    isMark(C)    (((C) == 'M' || (C) == 'm')? TRUE : FALSE)
  209. #define    isSpace(C)    (((C) == 'S' || (C) == 's')? TRUE : FALSE)
  210. #define    MarkLength(C)    (C) == 'M'? 1000 : 250
  211. #define    SpaceLength(C)    (C) == 'S'? 1000 : 250
  212. #define    MAX_MARKS    32
  213.  
  214. void
  215. DecodeLineStyle(StyleStr, Mark, Space, N)
  216. char    *StyleStr;
  217. int    **Mark;
  218. int    **Space;
  219. int    *N;
  220. {
  221.     register char    *S;
  222.     int    macc, sacc;
  223.     int    i, j, k;
  224.     int    mark[MAX_MARKS];
  225.     int    space[MAX_MARKS];
  226.     bool    ProcessingMark;
  227.  
  228.     /*
  229.      * 'S' = 1000   micron space
  230.      * 's' =  250   micron space
  231.      * 'M' = 1000     micron line
  232.      * 'm' =  250    micron line
  233.      */
  234.     ProcessingMark = TRUE; /* to satisfy some compilers */
  235.     if (isMark(*StyleStr))
  236.         ProcessingMark = TRUE;
  237.     else if (isSpace(*StyleStr))
  238.         ProcessingMark = FALSE;
  239.     else {
  240.         fprintf(stderr,"(DecodeLineStyle) Invalid Line style: \"%s\"\n", StyleStr);
  241.         ErrorExit();
  242.     }
  243.     memset((char *)mark, 0, MAX_MARKS * sizeof(int));
  244.     memset((char *)space, 0, MAX_MARKS * sizeof(int));
  245.     for (macc=0,sacc=0,i=0,j=0,S=StyleStr ; *S && i<MAX_MARKS && j<MAX_MARKS; ++S) {
  246.         if (isMark(*S) && ProcessingMark) {
  247.             macc += MarkLength(*S);
  248.         } else if (isMark(*S) && !ProcessingMark) {
  249.             space[j++] = sacc;
  250.             macc = MarkLength(*S);
  251.             ProcessingMark = TRUE;
  252.         } else if (isSpace(*S) && ProcessingMark) {
  253.             mark[i++] = macc;
  254.             sacc = SpaceLength(*S);
  255.             ProcessingMark = FALSE;
  256.         } else if (isSpace(*S) && !ProcessingMark) {
  257.             sacc += SpaceLength(*S);
  258.         } else if (isspace(*S)) {
  259.             continue;
  260.         } else {
  261.             fprintf(stderr,"(DecodeLineStyle) Invalid line style code '%c' in \"%s\"\n", *S, StyleStr);
  262.             ErrorExit();
  263.         }
  264.     }
  265.     if (ProcessingMark)
  266.         mark[i++] = macc;
  267.     else
  268.         space[j++] = sacc;
  269.  
  270.     k = max(i,j);
  271.  
  272.     if (((*Mark) = (int *)calloc(k, sizeof(int))) == (int *)NULL) {
  273.         perror("(DecodeLineStyle) ");
  274.         ErrorExit();
  275.     }
  276.     if (((*Space) = (int *)calloc(k, sizeof(int))) == (int *)NULL) {
  277.         perror("(DecodeLineStyle) ");
  278.         ErrorExit();
  279.     }
  280.     for (i=0; i<k; i++) {
  281.         (*Mark)[i] = mark[i];
  282.         (*Space)[i] = space[i];
  283.     }
  284.     *N = k;
  285. }
  286.  
  287. void
  288. InitLineStyles()
  289. {
  290.     int        i;
  291.     LATOM        *A;
  292.     LINESTYLE    *LS;
  293.     int        *Mark, *Space;
  294.     int        N;
  295. #ifdef    ANSI_C
  296.     LINESTYLE    *LSAlloc(int , int *, int *);
  297. #else
  298.     LINESTYLE    *LSAlloc();
  299. #endif
  300.  
  301.     LineStyleList.list_n = 0;
  302.     LineStyleList.list_head = (LATOM *)NULL;
  303.     LineStyleList.list_tail = (LATOM *)NULL;
  304.  
  305.     LS = LSAlloc(0, (int *)NULL, (int *)NULL);
  306.     Append(&LineStyleList, (generic *)LS);
  307.  
  308.     for (i=1,A=V_LineStyle.val_u.udt_set.list_head; A; A=A->la_next,i++) {
  309.         DecodeLineStyle(
  310.             (char *)A->la_p,
  311.             &Mark,
  312.             &Space,
  313.             &N
  314.         );
  315.         LS = LSAlloc(N, Mark, Space);
  316.         Append(&LineStyleList, (generic *)LS);
  317.     }
  318.  
  319.     NLineStyles = i;
  320.     if ((LineStyles = (LINESTYLE *)calloc(NLineStyles, sizeof(LINESTYLE))) == (LINESTYLE *)NULL) {
  321.         perror("(InitLineStyles) ");
  322.         ErrorExit();
  323.     }
  324.     for (i=0,A=LineStyleList.list_head; A; A=A->la_next,i++) {
  325.         LineStyles[i].ls_mark = ((LINESTYLE *)(A->la_p))->ls_mark;
  326.         LineStyles[i].ls_space = ((LINESTYLE *)(A->la_p))->ls_space;
  327.         LineStyles[i].ls_nelms = ((LINESTYLE *)(A->la_p))->ls_nelms;
  328.     }
  329. }
  330.  
  331.  
  332. LINESTYLE    *
  333. LSAlloc(n, Mark, Space)
  334. int    n;
  335. int    *Mark, *Space;
  336. {
  337.     register LINESTYLE    *LS;
  338.  
  339.     if ((LS = (LINESTYLE *)calloc(1, sizeof(LINESTYLE))) == (LINESTYLE *)NULL) {
  340.         perror("(LSAlloc) ");
  341.         assert(LS != (LINESTYLE *)NULL);
  342.         exit(errno);
  343.     }
  344.     LS->ls_mark = Mark;
  345.     LS->ls_space = Space;
  346.     LS->ls_nelms = n;
  347.     return(LS);
  348. }
  349.  
  350. void
  351. InitPointCodes()
  352. {
  353.     register LATOM    *A;
  354.     register int    i;
  355.     extern int    *PointCodes, NPointCodes;
  356.     extern VALUE    V_PointSymbol;
  357.  
  358.     if (!isSet(V_PointSymbol)) {
  359.         NPointCodes = 0;
  360.         return;
  361.     }
  362.  
  363.     /* count # codes    */
  364.     for (i=0,A=V_PointSymbol.val_u.udt_set.list_head; A; A=A->la_next, i++)
  365.         ;
  366.  
  367.     NPointCodes = i;
  368.     if (NPointCodes == 0)
  369.         return;
  370.  
  371.     if ((PointCodes = (int *)calloc(NPointCodes, sizeof(int))) == (int *)NULL) {
  372.         perror("(InitPointCodes) ");
  373.         assert(PointCodes != (int *)NULL);
  374.         ErrorExit();
  375.     }
  376.     for (i=0,A=V_PointSymbol.val_u.udt_set.list_head; A; A=A->la_next, i++) {
  377.         PointCodes[i] = atoi((char *)A->la_p);
  378.     }
  379. }
  380.  
  381.  
  382. void
  383. Plot2D(Data, NTuples, TupleSize)
  384. FLOAT    **Data;
  385. int    NTuples;
  386. int    TupleSize;
  387. {
  388.     int    i, j;
  389.     FLOAT    *X, *Y;
  390.     FLOAT    x, y;
  391.     FLOAT    x0, y0;
  392.     FLOAT    xtick, ytick;
  393.     FLOAT    Xmin,Xmax, Ymin,Ymax;
  394.     FLOAT    MedianX, StdDevX;
  395.     FLOAT    MedianY, StdDevY;
  396.     int    nxsub, nysub;
  397.     bool    AutoRange, AutoDomain;
  398.     bool    LogX, LogY;
  399.     char    Xopt[16], Yopt[16];
  400. #ifdef    ANSI_C
  401.     void    StatisticsOfX(
  402.             FLOAT **,
  403.             int,
  404.             int ,
  405.             FLOAT *,
  406.             FLOAT *
  407.         );
  408.     void    StatisticsOfY(
  409.             FLOAT **,
  410.             int,
  411.             int ,
  412.             FLOAT *,
  413.             FLOAT *
  414.         );
  415.     double    log10(double);
  416. #else
  417.     void    StatisticsOfX();
  418.     void    StatisticsOfY();
  419.     double    log10();
  420. #endif
  421.  
  422.     pladv(0);
  423.     ViewPort = VtoRect(V_ViewPort);
  424.  
  425.     if (isString(V_AspectRatio)) {
  426.         /* automatic --> STRETCH */
  427.         plvpor(ViewPort.XMin, ViewPort.XMax, ViewPort.YMin, ViewPort.YMax);
  428. /*debug
  429.         plvsta();
  430. debug*/
  431.     } else {
  432.         assert(isDbl(V_AspectRatio));
  433.         plvasp(VtoDbl(V_AspectRatio));
  434.     }
  435.     
  436.     LogX = (Vstrcmp(V_PlotType, "loglin") || Vstrcmp(V_PlotType,"loglog"))?
  437.             TRUE : FALSE;
  438.     LogY = (Vstrcmp(V_PlotType, "linlog") || Vstrcmp(V_PlotType,"loglog"))?
  439.             TRUE : FALSE;
  440.  
  441.     /* window bounds */
  442.     AutoDomain = FALSE;
  443.     if (isString(V_Domain) && Vstrcmp(V_Domain, "All")) {
  444.         /* all points    */
  445.         Xmin = UserRect.XMin;
  446.         Xmax = UserRect.XMax;
  447.     } else if (isString(V_Domain) && Vstrcmp(V_Domain,"Automatic")) {
  448.         /* automatic     */
  449.         AutoDomain = TRUE;
  450.         StatisticsOfX(Data, NTuples, TupleSize, &MedianX, &StdDevX);
  451.         x = MedianX - 3.0*StdDevX;
  452.         Xmin = UserRect.XMin < x? x : UserRect.XMin;
  453.         x = MedianX + 3.0*StdDevX;
  454.         Xmax = UserRect.XMax > x? x : UserRect.XMax;
  455.     } else {
  456.         assert(isInterval(V_Domain));
  457.         Xmin = V_Domain.val_u.udt_interval.int_lo;
  458.         Xmax = V_Domain.val_u.udt_interval.int_hi;
  459.     }
  460.     if (LogX) {
  461.         if (Xmin == 0.0 || Xmax == 0.0) {
  462.             /* error */
  463.             fprintf(stderr,"ListPlot: Trying to find log10(0) in data! Stop.\n");
  464.             ErrorExit();
  465.         }
  466.         Xmin = log10((double)Xmin);
  467.         Xmax = log10((double)Xmax);
  468.     }
  469.  
  470.     AutoRange = FALSE;
  471.     if (isString(V_Range) && Vstrcmp(V_Range, "All")) {
  472.         /* all */
  473.         Ymin = UserRect.YMin;
  474.         Ymax = UserRect.YMax;
  475.     } else if (isString(V_Range) && Vstrcmp(V_Range,"Automatic")) {
  476.         /* automatic     */
  477.         AutoRange = TRUE;
  478.         StatisticsOfY(Data, NTuples, TupleSize, &MedianY, &StdDevY);
  479.         y = MedianY - 3.0*StdDevY;
  480.         Ymin = UserRect.YMin < y? y : UserRect.YMin;
  481.         y = MedianY + 3.0*StdDevY;
  482.         Ymax = UserRect.YMax > y? y : UserRect.YMax;
  483.     } else {
  484.         assert(isInterval(V_Range));
  485.         Ymin = V_Range.val_u.udt_interval.int_lo;
  486.         Ymax = V_Range.val_u.udt_interval.int_hi;
  487.     }
  488.     if (LogY) {
  489.         if (Ymin == 0.0 || Ymax == 0.0) {
  490.             /* error */
  491.             fprintf(stderr,"ListPlot: Trying to find log10(0) in data! Stop.\n");
  492.             ErrorExit();
  493.         }
  494.         Ymin = log10((double)Ymin);
  495.         Ymax = log10((double)Ymax);
  496.     }
  497.     plwind(
  498.         Xmin,
  499.         Xmax,
  500.         Ymin,
  501.         Ymax
  502.     );
  503.     plschr(0., VtoDbl(V_AnnotationScale));
  504.  
  505.  
  506.     /* X and Y ticks */
  507.     if (isString(V_XTick)) {
  508.         /* automatic */
  509.         xtick = 0.;
  510.         nxsub = 0;
  511.     } else {
  512.         xtick = V_XTick.val_u.udt_interval.int_lo;
  513.         if (LogX) xtick = log10((double)xtick);
  514.         nxsub = V_XTick.val_u.udt_interval.int_hi;
  515.     }
  516.     if (isString(V_YTick)) {
  517.         /* automatic */
  518.         ytick = 0.;
  519.         nysub = 0;
  520.     } else {
  521.         ytick = V_YTick.val_u.udt_interval.int_lo;
  522.         if (LogY) ytick = log10((double)ytick);
  523.         nysub = V_YTick.val_u.udt_interval.int_hi;
  524.     }
  525.  
  526.     /* Axis options    */
  527.     Xopt[0] = (char)NULL;
  528.     Yopt[0] = (char)NULL;
  529.     if (isBoolean(V_Boxed) && VtoBoolean(V_Boxed)) {
  530.         strcat(Xopt, "bc");
  531.         strcat(Yopt, "bc");
  532.     } else if (isString(V_Boxed)) {
  533.         if (Vstrchr(V_Boxed, 'b'))
  534.             strcat(Xopt, "b");
  535.         if (Vstrchr(V_Boxed, 't'))
  536.             strcat(Xopt, "c");
  537.         if (Vstrchr(V_Boxed, 'r'))
  538.             strcat(Yopt, "c");
  539.         if (Vstrchr(V_Boxed, 'l'))
  540.             strcat(Yopt, "b");
  541.     } else {
  542.         ;
  543.     }
  544.     if (LogX)
  545.         strcat(Xopt, "l");
  546.     if (LogY)
  547.         strcat(Yopt, "l");
  548.     strcat(Xopt, "nst");
  549.     strcat(Yopt, "nstv");
  550.  
  551.     if (isString(V_Origin) && Vstrcmp(V_Origin, "Median")) {
  552.         if (AutoDomain == FALSE) {
  553.             StatisticsOfX(Data,NTuples,TupleSize,&MedianX,&StdDevX);
  554.         }
  555.         if (AutoRange == FALSE) {
  556.             StatisticsOfY(Data,NTuples,TupleSize,&MedianY,&StdDevY);
  557.         }
  558.         plaxes(
  559.             MedianX,MedianY,
  560.             Xopt, xtick, nxsub,
  561.             Yopt, ytick,nysub
  562.         );
  563.     } else if (isInterval(V_Origin)) {
  564.         x0 = V_Origin.val_u.udt_interval.int_lo;
  565.         y0 = V_Origin.val_u.udt_interval.int_hi;
  566.         plaxes(
  567.             x0,y0,
  568.             Xopt, xtick, nxsub,
  569.             Yopt, ytick,nysub
  570.         );
  571.     } else if (isString(V_Origin) && Vstrcmp(V_Origin, "Automatic")) {
  572.         plbox(Xopt, xtick, nxsub, Yopt, ytick,nysub);
  573.     } else {
  574.         plbox(Xopt, xtick, nxsub, Yopt, ytick,nysub);
  575.     }
  576.  
  577.     /* gridding    */
  578.     if (VtoBoolean(V_Gridding)) {
  579.         /* gridding on */
  580.         plstyl(1, &mark1, &space1);
  581.         plbox("g", xtick, nxsub, "g", ytick,nysub);
  582.         plstyl(0, &mark0, &space0);
  583.     }
  584.  
  585.     /* plot curves    */
  586.     if ((X = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
  587.         perror("(Plot2D) ");
  588.         ErrorExit();
  589.     }
  590.     if ((Y = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
  591.         perror("(Plot2D) ");
  592.         ErrorExit();
  593.     }
  594.     if (VtoBoolean(V_SupplyAbscissa)) {
  595.         for (i=0; i<NTuples; i++)
  596.             X[i] = LogX? log10((double)(i+1)) : (FLOAT)i;
  597.     } else {
  598.         for (i=0; i<NTuples; i++)
  599.             X[i] = LogX? log10((double)Data[i][0]) : Data[i][0];
  600.     }
  601.     for (i=1; i<TupleSize; i++) {
  602.         if (NLineStyles > 0) {
  603.             /* use user-supplied patterns */
  604.             plstyl(
  605.                 LineStyles[(i-1)%NLineStyles].ls_nelms,
  606.                 LineStyles[(i-1)%NLineStyles].ls_mark,
  607.                 LineStyles[(i-1)%NLineStyles].ls_space
  608.             );
  609.         } else {
  610.             pllsty((i-1)%8 + 1);
  611.         }
  612.         plcol(i);
  613.         for (j=0; j<NTuples; j++)
  614.             Y[j] = LogY? log10((double)Data[j][i]) : Data[j][i];
  615.  
  616.         if (VtoBoolean(V_PlotJoined))
  617.             plline(NTuples, X, Y);
  618.  
  619.         if (VtoBoolean(V_PlotPoints)) {
  620.             plschr(0., VtoDbl(V_PointScale));
  621.             pllsty(1);
  622.             if (NPointCodes > 0) {
  623.                 plpoin(NTuples, X, Y, PointCodes[(i-1)%NPointCodes]);
  624.             } else {
  625.                 plpoin(NTuples, X, Y, i);
  626.             }
  627.         }
  628.     }
  629.  
  630.  
  631.     /* restore line styles etc... */
  632.     pllsty(1);
  633.     plcol(1);
  634.     plschr(0., VtoDbl(V_LabelScale));
  635.     pllab(VtoString(V_XLabel), VtoString(V_YLabel), "");
  636.     plschr(0., VtoDbl(V_TitleScale));
  637.     pllab("","", VtoString(V_PlotTitle));
  638. }
  639.  
  640. #define    DTR(D)    (double)((D)*3.1415927/180.0)
  641.  
  642. void
  643. PolarPlot(Data, NTuples, TupleSize)
  644. FLOAT    **Data;
  645. int    NTuples;
  646. int    TupleSize;
  647. {
  648.     int    i, j;
  649.     FLOAT    *A, *R;
  650.     FLOAT    x, y;
  651.     FLOAT    xtick, ytick;
  652.     FLOAT    atick, rtick;
  653.     FLOAT    Xmin,Xmax, Ymin,Ymax;
  654.     FLOAT    r, rmax, ang, amax;
  655.     FLOAT    MedianX, StdDevX;
  656.     FLOAT    MedianY, StdDevY;
  657.     FLOAT    da;
  658.     FLOAT    x0,y0, xn,yn;
  659.     double    pi;
  660.     char    s[32];
  661.     int    nxsub, nysub, nrsub, nasub;
  662.     bool    InRadians, XisAngular;
  663.     bool    AutoDomain, AutoRange;
  664. #ifdef    ANSI_C
  665.     void    StatisticsOfX(
  666.             FLOAT **,
  667.             int ,
  668.             int ,
  669.             FLOAT *,
  670.             FLOAT *
  671.         );
  672.     void    StatisticsOfY(
  673.             FLOAT **,
  674.             int ,
  675.             int ,
  676.             FLOAT *,
  677.             FLOAT *
  678.         );
  679.     double    log10(double); double fabs();
  680.     double    sqrt(), sin(), cos(), atan();
  681. #else
  682.     void    StatisticsOfX();
  683.     void    StatisticsOfY();
  684.     double    log10(); double fabs();
  685.     double    sqrt(), sin(), cos(), atan();
  686. #endif
  687.  
  688.     pi = 4.0 * atan((double)1.0);
  689.     InRadians = Vstrcmp(V_AngularUnit, "degrees")? FALSE : TRUE;
  690.     XisAngular = Vstrcmp(V_PolarVariable,"angle")? TRUE : FALSE;
  691.  
  692.     pladv(0);
  693.  
  694.     if (isString(V_AspectRatio)) {
  695.         /* automatic --> STRETCH */
  696.         plvasp(1.0);
  697.     } else {
  698.         assert(isDbl(V_AspectRatio));
  699.         plvasp(VtoDbl(V_AspectRatio));
  700.     }
  701.  
  702.     /* window bounds */
  703.     AutoDomain = FALSE;
  704.     if (isString(V_Domain) && Vstrcmp(V_Domain, "All")) {
  705.         /* all points    */
  706.         Xmin = UserRect.XMin;
  707.         Xmax = UserRect.XMax;
  708.     } else if (isString(V_Domain) && Vstrcmp(V_Domain,"Automatic")) {
  709.         /* automatic     */
  710.         AutoDomain = TRUE;
  711.         StatisticsOfX(Data, NTuples, TupleSize, &MedianX, &StdDevX);
  712.         x = MedianX - 3.0*StdDevX;
  713.         Xmin = UserRect.XMin < x? x : UserRect.XMin;
  714.         x = MedianX + 3.0*StdDevX;
  715.         Xmax = UserRect.XMax > x? x : UserRect.XMax;
  716.     } else {
  717.         assert(isInterval(V_Domain));
  718.         Xmin = V_Domain.val_u.udt_interval.int_lo;
  719.         Xmax = V_Domain.val_u.udt_interval.int_hi;
  720.     }
  721.  
  722.     AutoRange = FALSE;
  723.     if (isString(V_Range) && Vstrcmp(V_Range, "All")) {
  724.         /* all */
  725.         Ymin = UserRect.YMin;
  726.         Ymax = UserRect.YMax;
  727.     } else if (isString(V_Range) && Vstrcmp(V_Range,"Automatic")) {
  728.         /* automatic     */
  729.         AutoRange = TRUE;
  730.         StatisticsOfY(Data, NTuples, TupleSize, &MedianY, &StdDevY);
  731.         y = MedianY - 3.0*StdDevY;
  732.         Ymin = UserRect.YMin < y? y : UserRect.YMin;
  733.         y = MedianY + 3.0*StdDevY;
  734.         Ymax = UserRect.YMax > y? y : UserRect.YMax;
  735.     } else {
  736.         assert(isInterval(V_Range));
  737.         Ymin = V_Range.val_u.udt_interval.int_lo;
  738.         Ymax = V_Range.val_u.udt_interval.int_hi;
  739.     }
  740.     rmax = XisAngular? Ymax : Xmax;
  741.     plwind(
  742.         -1.3*rmax,
  743.         1.3*rmax,
  744.         -1.3*rmax,
  745.         1.3*rmax
  746.     );
  747.  
  748.  
  749.     /* X and Y ticks */
  750.     if (isString(V_XTick)) {
  751.         /* automatic */
  752.         xtick = XisAngular? (InRadians? pi/6.0 : 30.0) : 0.25 * rmax;
  753.         nxsub = 4;
  754.     } else {
  755.         xtick = V_XTick.val_u.udt_interval.int_lo;
  756.         nxsub = V_XTick.val_u.udt_interval.int_hi;
  757.     }
  758.     if (isString(V_YTick)) {
  759.         /* automatic */
  760.         ytick = XisAngular? 0.25 * rmax : (InRadians? pi/6.0 : 30.0);
  761.         nysub = 4;
  762.     } else {
  763.         ytick = V_YTick.val_u.udt_interval.int_lo;
  764.         nysub = V_YTick.val_u.udt_interval.int_hi;
  765.     }
  766.     rtick = XisAngular? ytick : xtick;
  767.     atick = XisAngular? xtick : ytick;
  768.     nasub = XisAngular? nxsub : nysub;
  769.     nrsub = XisAngular? nysub : nxsub;
  770.     amax = InRadians? 2.0 * pi : 360.0;
  771.  
  772.     /* Axis options    */
  773.     if (VtoBoolean(V_Gridding)) {
  774.         /* draw circles    */
  775.         plstyl(0, &mark1, space1);
  776.         /*  circular ticks    */
  777.         for (i=0,r=(rtick/nrsub); r<rmax; r += (rtick/nrsub),i++) {
  778.             da = i%nrsub? 0.05 * atick : 0.1 * atick;
  779.             for (ang=0; ang<amax; ang += atick) {
  780.                 x0 = r * cos(InRadians? ang-da : DTR(ang-da)); 
  781.                 y0 = r * sin(InRadians? ang-da : DTR(ang-da)); 
  782.                 xn = r * cos(InRadians? ang+da : DTR(ang+da)); 
  783.                 yn = r * sin(InRadians? ang+da : DTR(ang+da)); 
  784.                 pljoin(x0,y0,xn,yn);
  785.             }
  786.         }
  787. #ifdef    CIRCULAR_TICKS
  788.         /* major circles    */
  789.         for (r=rtick; r<rmax; r += rtick) {
  790.             x0 = r; y0 = 0.0;
  791.             for (ang=1.0; ang<=360.0; ang++) {
  792.                 xn = r * cos((double)(DTR(ang))); 
  793.                 yn = r * sin((double)(DTR(ang))); 
  794.                 pljoin(x0,y0,xn,yn);
  795.                 x0 = xn; y0 = yn;
  796.             }
  797.         }
  798. #endif
  799.         /* draw spokes    */
  800.         for (i=0,ang=0.0; ang<amax; ang += (atick/nasub),i++) {
  801.             xn = rmax * cos(InRadians? ang : DTR(ang));
  802.             yn = rmax * sin(InRadians? ang : DTR(ang));
  803.             if (i%nasub == 0) {
  804.                 pljoin((FLOAT)0.0,(FLOAT)0.0,xn,yn);
  805.             } else {
  806.                 pljoin(xn,yn,1.05*xn,1.05*yn);
  807.             }
  808.         }
  809.     }
  810.     /* write angle labels    */
  811.     plschr(0., VtoDbl(V_AnnotationScale));
  812.     j = amax / atick;
  813.     for (ang=0.0,i=0; i<j; ang += atick,i++) {
  814.         if (InRadians) {
  815.             xn = rmax * cos((double)ang);
  816.             yn = rmax * sin((double)ang);
  817.             if (fabs(ang-(2.0*pi)) > 1.0e-2)
  818.                 sprintf(s, "%.2f #gp", (float)ang/pi);
  819.         } else {
  820.             xn = rmax * cos((double)DTR(ang));
  821.             yn = rmax * sin((double)DTR(ang));
  822.             sprintf(s, "%d", (int)(ang+0.5));
  823.         }
  824.         if (xn >= 0)
  825.             plptex(xn,yn,xn,yn,-0.15, s);
  826.         else
  827.             plptex(xn,yn,-xn,-yn,1.15, s);
  828.     }
  829.     plstyl(0, &mark0, &space0);
  830.  
  831.     /* plot curves    */
  832.     if ((A = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
  833.         perror("(Plot2D) ");
  834.         ErrorExit();
  835.     }
  836.     if ((R = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
  837.         perror("(Plot2D) ");
  838.         ErrorExit();
  839.     }
  840.     if (VtoBoolean(V_SupplyAbscissa)) {
  841.         for (i=0; i<NTuples; i++)
  842.             if (XisAngular)
  843.                 A[i] = i*2.0*pi/(NTuples-1);
  844.             else
  845.                 R[i] = i*rmax/(NTuples-1);
  846.     } else {
  847.         for (i=0; i<NTuples; i++)
  848.             if (XisAngular)
  849.                 A[i] = InRadians? Data[i][0] :  DTR(Data[i][0]);
  850.             else
  851.                 R[i] = Data[i][0];
  852.     }
  853.     for (i=1; i<TupleSize; i++) {
  854.         if (NLineStyles > 0) {
  855.             /* use user-supplied patterns */
  856.             plstyl(
  857.                 LineStyles[(i-1)%NLineStyles].ls_nelms,
  858.                 LineStyles[(i-1)%NLineStyles].ls_mark,
  859.                 LineStyles[(i-1)%NLineStyles].ls_space
  860.             );
  861.         } else {
  862.             pllsty((i-1)%8);
  863.         }
  864.         plcol(i);
  865.  
  866.         /* get and convert data if necessary    */
  867.         if (InRadians && !XisAngular)
  868.             for (j=0; j<NTuples; j++)
  869.                 A[j] = Data[j][i];
  870.         else if ((!InRadians) && !XisAngular)
  871.             for (j=0; j<NTuples; j++)
  872.                 A[j] = DTR(Data[j][i]);
  873.         else
  874.             for (j=0; j<NTuples; j++)
  875.                 R[j] = Data[j][i];
  876.  
  877.         /* convert to cartesian    A<->X  R<->Y*/
  878.         for (j=0; j<NTuples; j++) {
  879.             x0 = R[j] * cos(A[j]);
  880.             y0 = R[j] * sin(A[j]);
  881.             A[j] = x0;
  882.             R[j] = y0;
  883.         }
  884.         if (VtoBoolean(V_PlotJoined)) {
  885.             x0 = A[0];
  886.             y0 = R[0];
  887.             for (j=1; j<NTuples; j++) {
  888.                 pljoin(A[j-1], R[j-1], A[j],R[j]);
  889.             }
  890.         }
  891.  
  892.         if (VtoBoolean(V_PlotPoints)) {
  893.             plschr(0., VtoDbl(V_PointScale));
  894.             pllsty(1);
  895.             if (NPointCodes > 0) {
  896.                 plpoin(NTuples, A, R, PointCodes[(i-1)%NPointCodes]);
  897.             } else {
  898.                 plpoin(NTuples, A, R, i);
  899.             }
  900.         }
  901.     }
  902.  
  903.     /* restore line styles etc... */
  904.     pllsty(1);
  905.     plcol(1);
  906.     plschr(0., VtoDbl(V_LabelScale));
  907.     pllab(VtoString(V_XLabel), VtoString(V_YLabel), "");
  908.     plschr(0., VtoDbl(V_TitleScale));
  909.     pllab("","", VtoString(V_PlotTitle));
  910. }
  911.  
  912.  
  913. void
  914. StatisticsOfX(Data, NTuples,TupleSize, Median, StdDev)
  915. FLOAT    **Data;
  916. int    NTuples;
  917. int    TupleSize;
  918. FLOAT    *Median;
  919. FLOAT    *StdDev;
  920. {
  921. #ifdef    ANSI_C
  922.     FLOAT    MedianOfX();
  923.     FLOAT    StdDevOfX();
  924. /*debug
  925.     FLOAT    MedianOfX(FLOAT **Data,int NTuples,int TupleSize);
  926.     FLOAT    StdDevOfX(FLOAT **Data,int NTuples,int TupleSize);
  927. debug*/
  928. #else
  929.     FLOAT    MedianOfX();
  930.     FLOAT    StdDevOfX();
  931. #endif
  932.     *StdDev = StdDevOfX(Data,NTuples,TupleSize);    
  933.     *Median = MedianOfX(Data,NTuples,TupleSize);
  934. }
  935.  
  936.  
  937. void
  938. StatisticsOfY(Data, NTuples,TupleSize, Median, StdDev)
  939. FLOAT    **Data;
  940. int    NTuples;
  941. int    TupleSize;
  942. FLOAT    *Median;
  943. FLOAT    *StdDev;
  944. {
  945. #ifdef    ANSI_C
  946.     FLOAT    MedianOfY();
  947.     FLOAT    StdDevOfY();
  948. /*debug
  949.     FLOAT    MedianOfY(FLOAT **Data,int NTuples,int TupleSize);
  950.     FLOAT    StdDevOfY(FLOAT **Data,int NTuples,int TupleSize);
  951. debug*/
  952. #else
  953.     FLOAT    MedianOfY();
  954.     FLOAT    StdDevOfY();
  955. #endif
  956.     *StdDev = StdDevOfY(Data,NTuples,TupleSize);    
  957.     *Median = MedianOfY(Data,NTuples,TupleSize);
  958. }
  959.  
  960.  
  961. FLOAT
  962. StdDevOfX(Data,NTuples,TupleSize)
  963. FLOAT    **Data;
  964. int    NTuples;
  965. int    TupleSize;
  966. {
  967.     int    i;
  968.     FLOAT    StdDev, x;
  969.     FLOAT    var, sum, mean;
  970.     double sqrt(), fabs();
  971.  
  972.     if (VtoBoolean(V_SupplyAbscissa)) {
  973.         return((FLOAT)(NTuples*NTuples)/12.0);
  974.     }
  975.  
  976.     sum = 0.0;
  977.     for (i=0; i<NTuples; i++)
  978.         sum += Data[i][0];
  979.     mean = sum / NTuples;
  980.  
  981.     var = 0.0;
  982.     for (i=0; i<NTuples; i++) {
  983.         x = fabs((double)(Data[i][0] -  mean));
  984.         var += x * x;
  985.     }
  986.     var /= (NTuples-1);
  987.     StdDev = sqrt((double)var);
  988.     return(StdDev);    
  989. }
  990.  
  991.  
  992. FLOAT
  993. StdDevOfY(Data,NTuples,TupleSize)
  994. FLOAT    **Data;
  995. int    NTuples;
  996. int    TupleSize;
  997. {
  998.     int    i,j,D;
  999.     FLOAT    StdDev, x;
  1000.     FLOAT    var, sum, mean;
  1001.     FLOAT    n;
  1002.     double sqrt(), fabs();
  1003.  
  1004.     D = VtoBoolean(V_SupplyAbscissa) ? 0 : 1;
  1005.     n = NTuples * (TupleSize-D);
  1006.     sum = 0.0;
  1007.     for (i=0; i<NTuples; i++)
  1008.         for (j=D; j<TupleSize; j++) {
  1009.             sum += Data[i][j];
  1010.         }
  1011.     mean = sum / n;
  1012.  
  1013.     var = 0.0;
  1014.     for (i=0; i<NTuples; i++)
  1015.         for (j=D; j<TupleSize; j++) {
  1016.             x = fabs((double)(Data[i][j] -  mean));
  1017.             var += x * x;
  1018.         }
  1019.     var /= (n-1);
  1020.     StdDev = sqrt((double)var);
  1021.     return(StdDev);    
  1022. }
  1023.  
  1024. #define    BIG    1.0e30
  1025. #define    AFAC    1.5
  1026. #define    AMP    1.5
  1027.  
  1028. FLOAT
  1029. MedianOfX(Data,NTuples,TupleSize)
  1030. FLOAT    **Data;
  1031. int    NTuples;
  1032. int    TupleSize;
  1033. /* an iterative computation of the median...
  1034.  * Code taken from Numerical Recipes..
  1035.  */
  1036. {
  1037.     int    np, nm, i;
  1038.     FLOAT    xx,xp,xm,sumx,sum,eps,stemp,dum,ap,am,aa,a;
  1039.     FLOAT    Median;
  1040.     double    fabs();
  1041.  
  1042.     if (VtoBoolean(V_SupplyAbscissa))
  1043.         return((FLOAT)(0.5*NTuples));
  1044.  
  1045.     /* first estimate    */
  1046.     a = 0.5 * (Data[0][0] + Data[NTuples-1][0]);
  1047.     eps = fabs((double)(Data[NTuples-1][0] - Data[0][0]));
  1048.  
  1049.     am = -(ap=BIG);
  1050.     for (;;) {
  1051.         sum = sumx = 0.0;
  1052.         np = nm = 0;        /* # point above and below median */
  1053.         xm = -(xp = BIG);
  1054.  
  1055.         for (i=0; i<NTuples; i++) {
  1056.             xx = Data[i][0];
  1057.             if (xx != a) {
  1058.                 if (xx > a) {
  1059.                     ++np;
  1060.                     if (xx < xp) xp = xx;
  1061.                 } else if (xx < a) {
  1062.                     ++nm;
  1063.                     if (xx > xm) xm = xx;
  1064.                 }
  1065.                 sum += dum=1.0/(eps+fabs((double)(xx-a)));
  1066.                 sumx += xx * dum;
  1067.             }
  1068.         }
  1069.  
  1070.         stemp = (sumx / sum) - a;
  1071.         if ((np-nm) >= 2) {
  1072.             /* guess is too low make another pass */
  1073.             am = a;
  1074.             aa = stemp < 0.0? xp : xp + stemp*AMP;
  1075.             if (aa > ap) aa = 0.5*(a+ap);
  1076.             eps = AFAC * fabs((double)(aa-a));
  1077.             a = aa;
  1078.         } else if ((nm-np) >= 2) {
  1079.             /* guess is too high make another pass */
  1080.             ap = a;
  1081.             aa = stemp > 0.0? xm : xm+stemp*AMP;
  1082.             if (aa < am) aa = 0.5*(a+am);
  1083.             eps = AFAC*fabs((double)(aa-a));
  1084.             a = aa;
  1085.         } else {    /* got it    */
  1086.             if (NTuples%2 == 0) {
  1087.                 Median = 0.5*(np==nm? xp+xm : np>nm? a+xp : xm+a);
  1088.             } else {
  1089.                 Median = np == nm? a : np>nm? xp : xm;
  1090.             }
  1091.             return(Median);
  1092.         }
  1093.     }
  1094. }
  1095.  
  1096. FLOAT
  1097. MedianOfY(Data,N,M)
  1098. FLOAT    **Data;
  1099. int    N;
  1100. int    M;
  1101. /* an iterative computation of the median...
  1102.  * Code taken from Numerical Recipes..
  1103.  */
  1104. {
  1105.     int    n, np, nm, i,j,D;
  1106.     FLOAT    xx,xp,xm,sumx,sum,eps,stemp,dum,ap,am,aa,a;
  1107.     FLOAT    Median;
  1108.     double    fabs();
  1109.  
  1110.     D = VtoBoolean(V_SupplyAbscissa) ? 0 : 1;
  1111.     n = N * (M-D);
  1112.  
  1113.     /* first estimate    */
  1114.     a = 0.5 * (Data[0][D] + Data[N-1][M-1]);
  1115.     eps = fabs((double)(Data[N-1][M-1] - Data[0][D]));
  1116.  
  1117.     am = -(ap=BIG);
  1118.     for (;;) {
  1119.         sum = sumx = 0.0;
  1120.         np = nm = 0;        /* # point above and below median */
  1121.         xm = -(xp = BIG);
  1122.  
  1123.         for (i=0; i<N; i++) {
  1124.             for (j=D; j<M; j++) {
  1125.                 xx = Data[i][j];
  1126.                 if (xx != a) {
  1127.                     if (xx > a) {
  1128.                         ++np;
  1129.                         if (xx < xp) xp = xx;
  1130.                     } else if (xx < a) {
  1131.                         ++nm;
  1132.                         if (xx > xm) xm = xx;
  1133.                     }
  1134.                     sum += dum=1.0/(eps+fabs((double)(xx-a)));
  1135.                     sumx += xx * dum;
  1136.                 }
  1137.             }
  1138.         }
  1139.  
  1140.         stemp = (sumx / sum) - a;
  1141.         if ((np-nm) >= 2) {
  1142.             /* guess is too low make another pass */
  1143.             am = a;
  1144.             aa = stemp < 0.0? xp : xp + stemp*AMP;
  1145.             if (aa > ap) aa = 0.5*(a+ap);
  1146.             eps = AFAC * fabs((double)(aa-a));
  1147.             a = aa;
  1148.         } else if ((nm-np) >= 2) {
  1149.             /* guess is too high make another pass */
  1150.             ap = a;
  1151.             aa = stemp > 0.0? xm : xm+stemp*AMP;
  1152.             if (aa < am) aa = 0.5*(a+am);
  1153.             eps = AFAC*fabs((double)(aa-a));
  1154.             a = aa;
  1155.         } else {    /* got it    */
  1156.             if (n%2 == 0) {
  1157.                 Median = 0.5*(np==nm? xp+xm : np>nm? a+xp : xm+a);
  1158.             } else {
  1159.                 Median = np == nm? a : np>nm? xp : xm;
  1160.             }
  1161.             return(Median);
  1162.         }
  1163.     }
  1164. }
  1165.