amiga-request@abcfd20.larc.nasa.gov (Amiga Sources/Binaries Moderator) (08/23/90)
Submitted-by: Anthony M. Richardson <amr@dukee.egr.duke.edu>
Posting-number: Volume 90, Issue 244
Archive-name: applications/listplot/part03
#!/bin/sh
# This is a shell archive. Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file". To overwrite existing
# files, type "sh file -c". You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g.. If this archive is complete, you
# will see the following message at the end:
# "End of archive 3 (of 3)."
# Contents: Csrc/plotting.c
# Wrapped by tadguy@abcfd20 on Wed Aug 22 21:03:22 1990
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'Csrc/plotting.c' -a "${1}" != "-c" ; then
echo shar: Will not clobber existing file \"'Csrc/plotting.c'\"
else
echo shar: Extracting \"'Csrc/plotting.c'\" \(25452 characters\)
sed "s/^X//" >'Csrc/plotting.c' <<'END_OF_FILE'
X/*
X * plotting.c
X *
X * plotting routines.
X *
X */
X#include <stdio.h>
X#include <errno.h>
Xextern int errno;
X#include <assert.h>
X#include <string.h>
X#include <ctype.h>
X#include <signal.h>
X#include "datatypes.h"
X
Xvoid ErrorExit();
X
X#ifndef min
X# define min(A,B) A<B? A : B
X#endif
X#ifndef max
X# define max(A,B) A>B? A : B
X#endif
X
Xextern bool GraphicsInProgress;
X
X/* GoldenRatio = (1+Sqrt(5))/2 */
X#define GOLDENRATIO 1.6180339887499
X#define ASPECTRATIO 1.0/GOLDENRATIO
Xextern VALUE V_AngularUnit;
Xextern VALUE V_AnnotationScale;
Xextern VALUE V_AspectRatio;
Xextern VALUE V_Boxed;
Xextern VALUE V_Domain;
Xextern VALUE V_Gridding;
Xextern VALUE V_LabelScale;
Xextern VALUE V_LineColor;
Xextern VALUE V_LineStyle;
Xextern VALUE V_Orientation;
Xextern VALUE V_Origin;
Xextern VALUE V_PlotColor;
Xextern VALUE V_PlotDevice;
Xextern VALUE V_PlotJoined;
Xextern VALUE V_PlotPoints;
Xextern VALUE V_PlotTitle;
Xextern VALUE V_PlotType;
Xextern VALUE V_PointScale;
Xextern VALUE V_PointSymbol;
Xextern VALUE V_PolarVariable;
Xextern VALUE V_Range;
Xextern VALUE V_SubPages;
Xextern VALUE V_SupplyAbscissa;
Xextern VALUE V_TitleScale;
Xextern VALUE V_UseInputFile;
Xextern VALUE V_Verbose;
Xextern VALUE V_ViewPort;
Xextern VALUE V_XLabel;
Xextern VALUE V_XTick;
Xextern VALUE V_YLabel;
Xextern VALUE V_YTick;
X
Xchar *DeviceTypes[] = {
X "amiga",
X "printer",
X "iff",
X "hp",
X "aegis",
X "postscript",
X (char *)NULL
X};
X
Xtypedef struct line_style {
X int ls_type;
X int ls_nelms;
X int *ls_space;
X int *ls_mark;
X } LINESTYLE;
XLIST LineStyleList;
XLINESTYLE *LineStyles; int NLineStyles;
Xint *PointCodes, NPointCodes;
Xint space0 = 0, mark0 = 0, space1 = 1500, mark1 = 1500;
X
XRECT UserRect;
XRECT ViewPort = { 0.1,0.1, 0.9,0.9 };
XRECT WorldRect = {0.0,0.0, 1.0,1.0};
X
X
Xvoid
XListPlot(Data, NTuples, TupleSize)
XFLOAT **Data;
Xint NTuples, TupleSize;
X{
X void InitializeDevice();
X void InitializePlottingEnv();
X#ifdef ANSI_C
X void Plot2D(FLOAT **, int , int );
X void PolarPlot(FLOAT **, int , int );
X#else
X void Plot2D();
X void PolarPlot();
X#endif
X
X InitializeDevice();
X GraphicsInProgress = TRUE;
X
X InitializePlottingEnv(Data, NTuples, TupleSize);
X
X if (Vstrcmp(V_PlotType,"polar")) {
X PolarPlot(Data, NTuples, TupleSize);
X } else {
X Plot2D(Data, NTuples, TupleSize);
X }
X plend();
X}
X
X
Xvoid
XInitializeDevice()
X{
X register int i;
X register char *DeviceStr;
X
X
X plorient(Vstrcmp(V_Orientation, "portrait")? 1 : 0);
X plselect(stdout); /* write to stdout */
X
X
X DeviceStr = V_PlotDevice.val_u.udt_string;
X for (i=0; DeviceTypes[i]; i++)
X if (strcmp(DeviceStr, DeviceTypes[i]) == 0) {
X plbeg(
X i+1, /* device code */
X (int)V_SubPages.val_u.udt_interval.int_lo, /* nx */
X (int)V_SubPages.val_u.udt_interval.int_hi /* ny */
X );
X return;
X }
X fprintf(stderr,"(InitializeDevice) Unknown device type: %s\n", DeviceStr);
X ErrorExit();
X}
X
Xvoid
XInitializePlottingEnv(Data, NTuples, TupleSize)
XFLOAT **Data;
Xint NTuples;
Xint TupleSize;
X{
X int i, j;
X FLOAT *Ptr;
X void InitLineStyles();
X void InitPointCodes();
X
X /* viewport */
X ViewPort = VtoRect(V_ViewPort);
X
X /* user window */
X UserRect.XMin = MAX_FLOAT;
X UserRect.YMin = MAX_FLOAT;
X UserRect.XMax = -MAX_FLOAT;
X UserRect.YMax = -MAX_FLOAT;
X
X if (VtoBoolean(V_SupplyAbscissa)) {
X UserRect.XMin = 0.0;
X UserRect.XMax = (FLOAT)(NTuples-1);
X } else {
X for (i=0; i<NTuples; i++) {
X if (Data[i][0] < UserRect.XMin)
X UserRect.XMin = Data[i][0];
X if (Data[i][0] > UserRect.XMax)
X UserRect.XMax = Data[i][0];
X }
X }
X for (i=0; i<NTuples; i++)
X for (j=TupleSize-1,Ptr= &Data[i][1]; j>0; --j,Ptr++) {
X if (*Ptr < UserRect.YMin)
X UserRect.YMin = *Ptr;
X if (*Ptr > UserRect.YMax)
X UserRect.YMax = *Ptr;
X }
X
X
X /* font scaling */
X
X /* plot type */
X
X /* grid */
X
X /* line styles */
X InitLineStyles();
X
X
X /* line colors */
X
X /* point symbols */
X InitPointCodes();
X
X /* window type */
X
X /* X label */
X
X /* Y label */
X
X /* Plot title */
X
X /* aspect ratio */
X}
X
X#define isMark(C) (((C) == 'M' || (C) == 'm')? TRUE : FALSE)
X#define isSpace(C) (((C) == 'S' || (C) == 's')? TRUE : FALSE)
X#define MarkLength(C) (C) == 'M'? 1000 : 250
X#define SpaceLength(C) (C) == 'S'? 1000 : 250
X#define MAX_MARKS 32
X
Xvoid
XDecodeLineStyle(StyleStr, Mark, Space, N)
Xchar *StyleStr;
Xint **Mark;
Xint **Space;
Xint *N;
X{
X register char *S;
X int macc, sacc;
X int i, j, k;
X int mark[MAX_MARKS];
X int space[MAX_MARKS];
X bool ProcessingMark;
X
X /*
X * 'S' = 1000 micron space
X * 's' = 250 micron space
X * 'M' = 1000 micron line
X * 'm' = 250 micron line
X */
X ProcessingMark = TRUE; /* to satisfy some compilers */
X if (isMark(*StyleStr))
X ProcessingMark = TRUE;
X else if (isSpace(*StyleStr))
X ProcessingMark = FALSE;
X else {
X fprintf(stderr,"(DecodeLineStyle) Invalid Line style: \"%s\"\n", StyleStr);
X ErrorExit();
X }
X memset((char *)mark, 0, MAX_MARKS * sizeof(int));
X memset((char *)space, 0, MAX_MARKS * sizeof(int));
X for (macc=0,sacc=0,i=0,j=0,S=StyleStr ; *S && i<MAX_MARKS && j<MAX_MARKS; ++S) {
X if (isMark(*S) && ProcessingMark) {
X macc += MarkLength(*S);
X } else if (isMark(*S) && !ProcessingMark) {
X space[j++] = sacc;
X macc = MarkLength(*S);
X ProcessingMark = TRUE;
X } else if (isSpace(*S) && ProcessingMark) {
X mark[i++] = macc;
X sacc = SpaceLength(*S);
X ProcessingMark = FALSE;
X } else if (isSpace(*S) && !ProcessingMark) {
X sacc += SpaceLength(*S);
X } else if (isspace(*S)) {
X continue;
X } else {
X fprintf(stderr,"(DecodeLineStyle) Invalid line style code '%c' in \"%s\"\n", *S, StyleStr);
X ErrorExit();
X }
X }
X if (ProcessingMark)
X mark[i++] = macc;
X else
X space[j++] = sacc;
X
X k = max(i,j);
X
X if (((*Mark) = (int *)calloc(k, sizeof(int))) == (int *)NULL) {
X perror("(DecodeLineStyle) ");
X ErrorExit();
X }
X if (((*Space) = (int *)calloc(k, sizeof(int))) == (int *)NULL) {
X perror("(DecodeLineStyle) ");
X ErrorExit();
X }
X for (i=0; i<k; i++) {
X (*Mark)[i] = mark[i];
X (*Space)[i] = space[i];
X }
X *N = k;
X}
X
Xvoid
XInitLineStyles()
X{
X int i;
X LATOM *A;
X LINESTYLE *LS;
X int *Mark, *Space;
X int N;
X#ifdef ANSI_C
X LINESTYLE *LSAlloc(int , int *, int *);
X#else
X LINESTYLE *LSAlloc();
X#endif
X
X LineStyleList.list_n = 0;
X LineStyleList.list_head = (LATOM *)NULL;
X LineStyleList.list_tail = (LATOM *)NULL;
X
X LS = LSAlloc(0, (int *)NULL, (int *)NULL);
X Append(&LineStyleList, (generic *)LS);
X
X for (i=1,A=V_LineStyle.val_u.udt_set.list_head; A; A=A->la_next,i++) {
X DecodeLineStyle(
X (char *)A->la_p,
X &Mark,
X &Space,
X &N
X );
X LS = LSAlloc(N, Mark, Space);
X Append(&LineStyleList, (generic *)LS);
X }
X
X NLineStyles = i;
X if ((LineStyles = (LINESTYLE *)calloc(NLineStyles, sizeof(LINESTYLE))) == (LINESTYLE *)NULL) {
X perror("(InitLineStyles) ");
X ErrorExit();
X }
X for (i=0,A=LineStyleList.list_head; A; A=A->la_next,i++) {
X LineStyles[i].ls_mark = ((LINESTYLE *)(A->la_p))->ls_mark;
X LineStyles[i].ls_space = ((LINESTYLE *)(A->la_p))->ls_space;
X LineStyles[i].ls_nelms = ((LINESTYLE *)(A->la_p))->ls_nelms;
X }
X}
X
X
XLINESTYLE *
XLSAlloc(n, Mark, Space)
Xint n;
Xint *Mark, *Space;
X{
X register LINESTYLE *LS;
X
X if ((LS = (LINESTYLE *)calloc(1, sizeof(LINESTYLE))) == (LINESTYLE *)NULL) {
X perror("(LSAlloc) ");
X assert(LS != (LINESTYLE *)NULL);
X exit(errno);
X }
X LS->ls_mark = Mark;
X LS->ls_space = Space;
X LS->ls_nelms = n;
X return(LS);
X}
X
Xvoid
XInitPointCodes()
X{
X register LATOM *A;
X register int i;
X extern int *PointCodes, NPointCodes;
X extern VALUE V_PointSymbol;
X
X if (!isSet(V_PointSymbol)) {
X NPointCodes = 0;
X return;
X }
X
X /* count # codes */
X for (i=0,A=V_PointSymbol.val_u.udt_set.list_head; A; A=A->la_next, i++)
X ;
X
X NPointCodes = i;
X if (NPointCodes == 0)
X return;
X
X if ((PointCodes = (int *)calloc(NPointCodes, sizeof(int))) == (int *)NULL) {
X perror("(InitPointCodes) ");
X assert(PointCodes != (int *)NULL);
X ErrorExit();
X }
X for (i=0,A=V_PointSymbol.val_u.udt_set.list_head; A; A=A->la_next, i++) {
X PointCodes[i] = atoi((char *)A->la_p);
X }
X}
X
X
Xvoid
XPlot2D(Data, NTuples, TupleSize)
XFLOAT **Data;
Xint NTuples;
Xint TupleSize;
X{
X int i, j;
X FLOAT *X, *Y;
X FLOAT x, y;
X FLOAT x0, y0;
X FLOAT xtick, ytick;
X FLOAT Xmin,Xmax, Ymin,Ymax;
X FLOAT MedianX, StdDevX;
X FLOAT MedianY, StdDevY;
X int nxsub, nysub;
X bool AutoRange, AutoDomain;
X bool LogX, LogY;
X char Xopt[16], Yopt[16];
X#ifdef ANSI_C
X void StatisticsOfX(
X FLOAT **,
X int,
X int ,
X FLOAT *,
X FLOAT *
X );
X void StatisticsOfY(
X FLOAT **,
X int,
X int ,
X FLOAT *,
X FLOAT *
X );
X double log10(double);
X#else
X void StatisticsOfX();
X void StatisticsOfY();
X double log10();
X#endif
X
X pladv(0);
X ViewPort = VtoRect(V_ViewPort);
X
X if (isString(V_AspectRatio)) {
X /* automatic --> STRETCH */
X plvpor(ViewPort.XMin, ViewPort.XMax, ViewPort.YMin, ViewPort.YMax);
X/*debug
X plvsta();
Xdebug*/
X } else {
X assert(isDbl(V_AspectRatio));
X plvasp(VtoDbl(V_AspectRatio));
X }
X
X LogX = (Vstrcmp(V_PlotType, "loglin") || Vstrcmp(V_PlotType,"loglog"))?
X TRUE : FALSE;
X LogY = (Vstrcmp(V_PlotType, "linlog") || Vstrcmp(V_PlotType,"loglog"))?
X TRUE : FALSE;
X
X /* window bounds */
X AutoDomain = FALSE;
X if (isString(V_Domain) && Vstrcmp(V_Domain, "All")) {
X /* all points */
X Xmin = UserRect.XMin;
X Xmax = UserRect.XMax;
X } else if (isString(V_Domain) && Vstrcmp(V_Domain,"Automatic")) {
X /* automatic */
X AutoDomain = TRUE;
X StatisticsOfX(Data, NTuples, TupleSize, &MedianX, &StdDevX);
X x = MedianX - 3.0*StdDevX;
X Xmin = UserRect.XMin < x? x : UserRect.XMin;
X x = MedianX + 3.0*StdDevX;
X Xmax = UserRect.XMax > x? x : UserRect.XMax;
X } else {
X assert(isInterval(V_Domain));
X Xmin = V_Domain.val_u.udt_interval.int_lo;
X Xmax = V_Domain.val_u.udt_interval.int_hi;
X }
X if (LogX) {
X if (Xmin == 0.0 || Xmax == 0.0) {
X /* error */
X fprintf(stderr,"ListPlot: Trying to find log10(0) in data! Stop.\n");
X ErrorExit();
X }
X Xmin = log10((double)Xmin);
X Xmax = log10((double)Xmax);
X }
X
X AutoRange = FALSE;
X if (isString(V_Range) && Vstrcmp(V_Range, "All")) {
X /* all */
X Ymin = UserRect.YMin;
X Ymax = UserRect.YMax;
X } else if (isString(V_Range) && Vstrcmp(V_Range,"Automatic")) {
X /* automatic */
X AutoRange = TRUE;
X StatisticsOfY(Data, NTuples, TupleSize, &MedianY, &StdDevY);
X y = MedianY - 3.0*StdDevY;
X Ymin = UserRect.YMin < y? y : UserRect.YMin;
X y = MedianY + 3.0*StdDevY;
X Ymax = UserRect.YMax > y? y : UserRect.YMax;
X } else {
X assert(isInterval(V_Range));
X Ymin = V_Range.val_u.udt_interval.int_lo;
X Ymax = V_Range.val_u.udt_interval.int_hi;
X }
X if (LogY) {
X if (Ymin == 0.0 || Ymax == 0.0) {
X /* error */
X fprintf(stderr,"ListPlot: Trying to find log10(0) in data! Stop.\n");
X ErrorExit();
X }
X Ymin = log10((double)Ymin);
X Ymax = log10((double)Ymax);
X }
X plwind(
X Xmin,
X Xmax,
X Ymin,
X Ymax
X );
X plschr(0., VtoDbl(V_AnnotationScale));
X
X
X /* X and Y ticks */
X if (isString(V_XTick)) {
X /* automatic */
X xtick = 0.;
X nxsub = 0;
X } else {
X xtick = V_XTick.val_u.udt_interval.int_lo;
X if (LogX) xtick = log10((double)xtick);
X nxsub = V_XTick.val_u.udt_interval.int_hi;
X }
X if (isString(V_YTick)) {
X /* automatic */
X ytick = 0.;
X nysub = 0;
X } else {
X ytick = V_YTick.val_u.udt_interval.int_lo;
X if (LogY) ytick = log10((double)ytick);
X nysub = V_YTick.val_u.udt_interval.int_hi;
X }
X
X /* Axis options */
X Xopt[0] = (char)NULL;
X Yopt[0] = (char)NULL;
X if (isBoolean(V_Boxed) && VtoBoolean(V_Boxed)) {
X strcat(Xopt, "bc");
X strcat(Yopt, "bc");
X } else if (isString(V_Boxed)) {
X if (Vstrchr(V_Boxed, 'b'))
X strcat(Xopt, "b");
X if (Vstrchr(V_Boxed, 't'))
X strcat(Xopt, "c");
X if (Vstrchr(V_Boxed, 'r'))
X strcat(Yopt, "c");
X if (Vstrchr(V_Boxed, 'l'))
X strcat(Yopt, "b");
X } else {
X ;
X }
X if (LogX)
X strcat(Xopt, "l");
X if (LogY)
X strcat(Yopt, "l");
X strcat(Xopt, "nst");
X strcat(Yopt, "nstv");
X
X if (isString(V_Origin) && Vstrcmp(V_Origin, "Median")) {
X if (AutoDomain == FALSE) {
X StatisticsOfX(Data,NTuples,TupleSize,&MedianX,&StdDevX);
X }
X if (AutoRange == FALSE) {
X StatisticsOfY(Data,NTuples,TupleSize,&MedianY,&StdDevY);
X }
X plaxes(
X MedianX,MedianY,
X Xopt, xtick, nxsub,
X Yopt, ytick,nysub
X );
X } else if (isInterval(V_Origin)) {
X x0 = V_Origin.val_u.udt_interval.int_lo;
X y0 = V_Origin.val_u.udt_interval.int_hi;
X plaxes(
X x0,y0,
X Xopt, xtick, nxsub,
X Yopt, ytick,nysub
X );
X } else if (isString(V_Origin) && Vstrcmp(V_Origin, "Automatic")) {
X plbox(Xopt, xtick, nxsub, Yopt, ytick,nysub);
X } else {
X plbox(Xopt, xtick, nxsub, Yopt, ytick,nysub);
X }
X
X /* gridding */
X if (VtoBoolean(V_Gridding)) {
X /* gridding on */
X plstyl(1, &mark1, &space1);
X plbox("g", xtick, nxsub, "g", ytick,nysub);
X plstyl(0, &mark0, &space0);
X }
X
X /* plot curves */
X if ((X = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
X perror("(Plot2D) ");
X ErrorExit();
X }
X if ((Y = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
X perror("(Plot2D) ");
X ErrorExit();
X }
X if (VtoBoolean(V_SupplyAbscissa)) {
X for (i=0; i<NTuples; i++)
X X[i] = LogX? log10((double)(i+1)) : (FLOAT)i;
X } else {
X for (i=0; i<NTuples; i++)
X X[i] = LogX? log10((double)Data[i][0]) : Data[i][0];
X }
X for (i=1; i<TupleSize; i++) {
X if (NLineStyles > 0) {
X /* use user-supplied patterns */
X plstyl(
X LineStyles[(i-1)%NLineStyles].ls_nelms,
X LineStyles[(i-1)%NLineStyles].ls_mark,
X LineStyles[(i-1)%NLineStyles].ls_space
X );
X } else {
X pllsty((i-1)%8 + 1);
X }
X plcol(i);
X for (j=0; j<NTuples; j++)
X Y[j] = LogY? log10((double)Data[j][i]) : Data[j][i];
X
X if (VtoBoolean(V_PlotJoined))
X plline(NTuples, X, Y);
X
X if (VtoBoolean(V_PlotPoints)) {
X plschr(0., VtoDbl(V_PointScale));
X pllsty(1);
X if (NPointCodes > 0) {
X plpoin(NTuples, X, Y, PointCodes[(i-1)%NPointCodes]);
X } else {
X plpoin(NTuples, X, Y, i);
X }
X }
X }
X
X
X /* restore line styles etc... */
X pllsty(1);
X plcol(1);
X plschr(0., VtoDbl(V_LabelScale));
X pllab(VtoString(V_XLabel), VtoString(V_YLabel), "");
X plschr(0., VtoDbl(V_TitleScale));
X pllab("","", VtoString(V_PlotTitle));
X}
X
X#define DTR(D) (double)((D)*3.1415927/180.0)
X
Xvoid
XPolarPlot(Data, NTuples, TupleSize)
XFLOAT **Data;
Xint NTuples;
Xint TupleSize;
X{
X int i, j;
X FLOAT *A, *R;
X FLOAT x, y;
X FLOAT xtick, ytick;
X FLOAT atick, rtick;
X FLOAT Xmin,Xmax, Ymin,Ymax;
X FLOAT r, rmax, ang, amax;
X FLOAT MedianX, StdDevX;
X FLOAT MedianY, StdDevY;
X FLOAT da;
X FLOAT x0,y0, xn,yn;
X double pi;
X char s[32];
X int nxsub, nysub, nrsub, nasub;
X bool InRadians, XisAngular;
X bool AutoDomain, AutoRange;
X#ifdef ANSI_C
X void StatisticsOfX(
X FLOAT **,
X int ,
X int ,
X FLOAT *,
X FLOAT *
X );
X void StatisticsOfY(
X FLOAT **,
X int ,
X int ,
X FLOAT *,
X FLOAT *
X );
X double log10(double); double fabs();
X double sqrt(), sin(), cos(), atan();
X#else
X void StatisticsOfX();
X void StatisticsOfY();
X double log10(); double fabs();
X double sqrt(), sin(), cos(), atan();
X#endif
X
X pi = 4.0 * atan((double)1.0);
X InRadians = Vstrcmp(V_AngularUnit, "degrees")? FALSE : TRUE;
X XisAngular = Vstrcmp(V_PolarVariable,"angle")? TRUE : FALSE;
X
X pladv(0);
X
X if (isString(V_AspectRatio)) {
X /* automatic --> STRETCH */
X plvasp(1.0);
X } else {
X assert(isDbl(V_AspectRatio));
X plvasp(VtoDbl(V_AspectRatio));
X }
X
X /* window bounds */
X AutoDomain = FALSE;
X if (isString(V_Domain) && Vstrcmp(V_Domain, "All")) {
X /* all points */
X Xmin = UserRect.XMin;
X Xmax = UserRect.XMax;
X } else if (isString(V_Domain) && Vstrcmp(V_Domain,"Automatic")) {
X /* automatic */
X AutoDomain = TRUE;
X StatisticsOfX(Data, NTuples, TupleSize, &MedianX, &StdDevX);
X x = MedianX - 3.0*StdDevX;
X Xmin = UserRect.XMin < x? x : UserRect.XMin;
X x = MedianX + 3.0*StdDevX;
X Xmax = UserRect.XMax > x? x : UserRect.XMax;
X } else {
X assert(isInterval(V_Domain));
X Xmin = V_Domain.val_u.udt_interval.int_lo;
X Xmax = V_Domain.val_u.udt_interval.int_hi;
X }
X
X AutoRange = FALSE;
X if (isString(V_Range) && Vstrcmp(V_Range, "All")) {
X /* all */
X Ymin = UserRect.YMin;
X Ymax = UserRect.YMax;
X } else if (isString(V_Range) && Vstrcmp(V_Range,"Automatic")) {
X /* automatic */
X AutoRange = TRUE;
X StatisticsOfY(Data, NTuples, TupleSize, &MedianY, &StdDevY);
X y = MedianY - 3.0*StdDevY;
X Ymin = UserRect.YMin < y? y : UserRect.YMin;
X y = MedianY + 3.0*StdDevY;
X Ymax = UserRect.YMax > y? y : UserRect.YMax;
X } else {
X assert(isInterval(V_Range));
X Ymin = V_Range.val_u.udt_interval.int_lo;
X Ymax = V_Range.val_u.udt_interval.int_hi;
X }
X rmax = XisAngular? Ymax : Xmax;
X plwind(
X -1.3*rmax,
X 1.3*rmax,
X -1.3*rmax,
X 1.3*rmax
X );
X
X
X /* X and Y ticks */
X if (isString(V_XTick)) {
X /* automatic */
X xtick = XisAngular? (InRadians? pi/6.0 : 30.0) : 0.25 * rmax;
X nxsub = 4;
X } else {
X xtick = V_XTick.val_u.udt_interval.int_lo;
X nxsub = V_XTick.val_u.udt_interval.int_hi;
X }
X if (isString(V_YTick)) {
X /* automatic */
X ytick = XisAngular? 0.25 * rmax : (InRadians? pi/6.0 : 30.0);
X nysub = 4;
X } else {
X ytick = V_YTick.val_u.udt_interval.int_lo;
X nysub = V_YTick.val_u.udt_interval.int_hi;
X }
X rtick = XisAngular? ytick : xtick;
X atick = XisAngular? xtick : ytick;
X nasub = XisAngular? nxsub : nysub;
X nrsub = XisAngular? nysub : nxsub;
X amax = InRadians? 2.0 * pi : 360.0;
X
X /* Axis options */
X if (VtoBoolean(V_Gridding)) {
X /* draw circles */
X plstyl(0, &mark1, space1);
X /* circular ticks */
X for (i=0,r=(rtick/nrsub); r<rmax; r += (rtick/nrsub),i++) {
X da = i%nrsub? 0.05 * atick : 0.1 * atick;
X for (ang=0; ang<amax; ang += atick) {
X x0 = r * cos(InRadians? ang-da : DTR(ang-da));
X y0 = r * sin(InRadians? ang-da : DTR(ang-da));
X xn = r * cos(InRadians? ang+da : DTR(ang+da));
X yn = r * sin(InRadians? ang+da : DTR(ang+da));
X pljoin(x0,y0,xn,yn);
X }
X }
X#ifdef CIRCULAR_TICKS
X /* major circles */
X for (r=rtick; r<rmax; r += rtick) {
X x0 = r; y0 = 0.0;
X for (ang=1.0; ang<=360.0; ang++) {
X xn = r * cos((double)(DTR(ang)));
X yn = r * sin((double)(DTR(ang)));
X pljoin(x0,y0,xn,yn);
X x0 = xn; y0 = yn;
X }
X }
X#endif
X /* draw spokes */
X for (i=0,ang=0.0; ang<amax; ang += (atick/nasub),i++) {
X xn = rmax * cos(InRadians? ang : DTR(ang));
X yn = rmax * sin(InRadians? ang : DTR(ang));
X if (i%nasub == 0) {
X pljoin((FLOAT)0.0,(FLOAT)0.0,xn,yn);
X } else {
X pljoin(xn,yn,1.05*xn,1.05*yn);
X }
X }
X }
X /* write angle labels */
X plschr(0., VtoDbl(V_AnnotationScale));
X j = amax / atick;
X for (ang=0.0,i=0; i<j; ang += atick,i++) {
X if (InRadians) {
X xn = rmax * cos((double)ang);
X yn = rmax * sin((double)ang);
X if (fabs(ang-(2.0*pi)) > 1.0e-2)
X sprintf(s, "%.2f #gp", (float)ang/pi);
X } else {
X xn = rmax * cos((double)DTR(ang));
X yn = rmax * sin((double)DTR(ang));
X sprintf(s, "%d", (int)(ang+0.5));
X }
X if (xn >= 0)
X plptex(xn,yn,xn,yn,-0.15, s);
X else
X plptex(xn,yn,-xn,-yn,1.15, s);
X }
X plstyl(0, &mark0, &space0);
X
X /* plot curves */
X if ((A = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
X perror("(Plot2D) ");
X ErrorExit();
X }
X if ((R = (FLOAT *)calloc(NTuples, sizeof(FLOAT))) == (FLOAT *)NULL) {
X perror("(Plot2D) ");
X ErrorExit();
X }
X if (VtoBoolean(V_SupplyAbscissa)) {
X for (i=0; i<NTuples; i++)
X if (XisAngular)
X A[i] = i*2.0*pi/(NTuples-1);
X else
X R[i] = i*rmax/(NTuples-1);
X } else {
X for (i=0; i<NTuples; i++)
X if (XisAngular)
X A[i] = InRadians? Data[i][0] : DTR(Data[i][0]);
X else
X R[i] = Data[i][0];
X }
X for (i=1; i<TupleSize; i++) {
X if (NLineStyles > 0) {
X /* use user-supplied patterns */
X plstyl(
X LineStyles[(i-1)%NLineStyles].ls_nelms,
X LineStyles[(i-1)%NLineStyles].ls_mark,
X LineStyles[(i-1)%NLineStyles].ls_space
X );
X } else {
X pllsty((i-1)%8);
X }
X plcol(i);
X
X /* get and convert data if necessary */
X if (InRadians && !XisAngular)
X for (j=0; j<NTuples; j++)
X A[j] = Data[j][i];
X else if ((!InRadians) && !XisAngular)
X for (j=0; j<NTuples; j++)
X A[j] = DTR(Data[j][i]);
X else
X for (j=0; j<NTuples; j++)
X R[j] = Data[j][i];
X
X /* convert to cartesian A<->X R<->Y*/
X for (j=0; j<NTuples; j++) {
X x0 = R[j] * cos(A[j]);
X y0 = R[j] * sin(A[j]);
X A[j] = x0;
X R[j] = y0;
X }
X if (VtoBoolean(V_PlotJoined)) {
X x0 = A[0];
X y0 = R[0];
X for (j=1; j<NTuples; j++) {
X pljoin(A[j-1], R[j-1], A[j],R[j]);
X }
X }
X
X if (VtoBoolean(V_PlotPoints)) {
X plschr(0., VtoDbl(V_PointScale));
X pllsty(1);
X if (NPointCodes > 0) {
X plpoin(NTuples, A, R, PointCodes[(i-1)%NPointCodes]);
X } else {
X plpoin(NTuples, A, R, i);
X }
X }
X }
X
X /* restore line styles etc... */
X pllsty(1);
X plcol(1);
X plschr(0., VtoDbl(V_LabelScale));
X pllab(VtoString(V_XLabel), VtoString(V_YLabel), "");
X plschr(0., VtoDbl(V_TitleScale));
X pllab("","", VtoString(V_PlotTitle));
X}
X
X
Xvoid
XStatisticsOfX(Data, NTuples,TupleSize, Median, StdDev)
XFLOAT **Data;
Xint NTuples;
Xint TupleSize;
XFLOAT *Median;
XFLOAT *StdDev;
X{
X#ifdef ANSI_C
X FLOAT MedianOfX();
X FLOAT StdDevOfX();
X/*debug
X FLOAT MedianOfX(FLOAT **Data,int NTuples,int TupleSize);
X FLOAT StdDevOfX(FLOAT **Data,int NTuples,int TupleSize);
Xdebug*/
X#else
X FLOAT MedianOfX();
X FLOAT StdDevOfX();
X#endif
X *StdDev = StdDevOfX(Data,NTuples,TupleSize);
X *Median = MedianOfX(Data,NTuples,TupleSize);
X}
X
X
Xvoid
XStatisticsOfY(Data, NTuples,TupleSize, Median, StdDev)
XFLOAT **Data;
Xint NTuples;
Xint TupleSize;
XFLOAT *Median;
XFLOAT *StdDev;
X{
X#ifdef ANSI_C
X FLOAT MedianOfY();
X FLOAT StdDevOfY();
X/*debug
X FLOAT MedianOfY(FLOAT **Data,int NTuples,int TupleSize);
X FLOAT StdDevOfY(FLOAT **Data,int NTuples,int TupleSize);
Xdebug*/
X#else
X FLOAT MedianOfY();
X FLOAT StdDevOfY();
X#endif
X *StdDev = StdDevOfY(Data,NTuples,TupleSize);
X *Median = MedianOfY(Data,NTuples,TupleSize);
X}
X
X
XFLOAT
XStdDevOfX(Data,NTuples,TupleSize)
XFLOAT **Data;
Xint NTuples;
Xint TupleSize;
X{
X int i;
X FLOAT StdDev, x;
X FLOAT var, sum, mean;
X double sqrt(), fabs();
X
X if (VtoBoolean(V_SupplyAbscissa)) {
X return((FLOAT)(NTuples*NTuples)/12.0);
X }
X
X sum = 0.0;
X for (i=0; i<NTuples; i++)
X sum += Data[i][0];
X mean = sum / NTuples;
X
X var = 0.0;
X for (i=0; i<NTuples; i++) {
X x = fabs((double)(Data[i][0] - mean));
X var += x * x;
X }
X var /= (NTuples-1);
X StdDev = sqrt((double)var);
X return(StdDev);
X}
X
X
XFLOAT
XStdDevOfY(Data,NTuples,TupleSize)
XFLOAT **Data;
Xint NTuples;
Xint TupleSize;
X{
X int i,j,D;
X FLOAT StdDev, x;
X FLOAT var, sum, mean;
X FLOAT n;
X double sqrt(), fabs();
X
X D = VtoBoolean(V_SupplyAbscissa) ? 0 : 1;
X n = NTuples * (TupleSize-D);
X sum = 0.0;
X for (i=0; i<NTuples; i++)
X for (j=D; j<TupleSize; j++) {
X sum += Data[i][j];
X }
X mean = sum / n;
X
X var = 0.0;
X for (i=0; i<NTuples; i++)
X for (j=D; j<TupleSize; j++) {
X x = fabs((double)(Data[i][j] - mean));
X var += x * x;
X }
X var /= (n-1);
X StdDev = sqrt((double)var);
X return(StdDev);
X}
X
X#define BIG 1.0e30
X#define AFAC 1.5
X#define AMP 1.5
X
XFLOAT
XMedianOfX(Data,NTuples,TupleSize)
XFLOAT **Data;
Xint NTuples;
Xint TupleSize;
X/* an iterative computation of the median...
X * Code taken from Numerical Recipes..
X */
X{
X int np, nm, i;
X FLOAT xx,xp,xm,sumx,sum,eps,stemp,dum,ap,am,aa,a;
X FLOAT Median;
X double fabs();
X
X if (VtoBoolean(V_SupplyAbscissa))
X return((FLOAT)(0.5*NTuples));
X
X /* first estimate */
X a = 0.5 * (Data[0][0] + Data[NTuples-1][0]);
X eps = fabs((double)(Data[NTuples-1][0] - Data[0][0]));
X
X am = -(ap=BIG);
X for (;;) {
X sum = sumx = 0.0;
X np = nm = 0; /* # point above and below median */
X xm = -(xp = BIG);
X
X for (i=0; i<NTuples; i++) {
X xx = Data[i][0];
X if (xx != a) {
X if (xx > a) {
X ++np;
X if (xx < xp) xp = xx;
X } else if (xx < a) {
X ++nm;
X if (xx > xm) xm = xx;
X }
X sum += dum=1.0/(eps+fabs((double)(xx-a)));
X sumx += xx * dum;
X }
X }
X
X stemp = (sumx / sum) - a;
X if ((np-nm) >= 2) {
X /* guess is too low make another pass */
X am = a;
X aa = stemp < 0.0? xp : xp + stemp*AMP;
X if (aa > ap) aa = 0.5*(a+ap);
X eps = AFAC * fabs((double)(aa-a));
X a = aa;
X } else if ((nm-np) >= 2) {
X /* guess is too high make another pass */
X ap = a;
X aa = stemp > 0.0? xm : xm+stemp*AMP;
X if (aa < am) aa = 0.5*(a+am);
X eps = AFAC*fabs((double)(aa-a));
X a = aa;
X } else { /* got it */
X if (NTuples%2 == 0) {
X Median = 0.5*(np==nm? xp+xm : np>nm? a+xp : xm+a);
X } else {
X Median = np == nm? a : np>nm? xp : xm;
X }
X return(Median);
X }
X }
X}
X
XFLOAT
XMedianOfY(Data,N,M)
XFLOAT **Data;
Xint N;
Xint M;
X/* an iterative computation of the median...
X * Code taken from Numerical Recipes..
X */
X{
X int n, np, nm, i,j,D;
X FLOAT xx,xp,xm,sumx,sum,eps,stemp,dum,ap,am,aa,a;
X FLOAT Median;
X double fabs();
X
X D = VtoBoolean(V_SupplyAbscissa) ? 0 : 1;
X n = N * (M-D);
X
X /* first estimate */
X a = 0.5 * (Data[0][D] + Data[N-1][M-1]);
X eps = fabs((double)(Data[N-1][M-1] - Data[0][D]));
X
X am = -(ap=BIG);
X for (;;) {
X sum = sumx = 0.0;
X np = nm = 0; /* # point above and below median */
X xm = -(xp = BIG);
X
X for (i=0; i<N; i++) {
X for (j=D; j<M; j++) {
X xx = Data[i][j];
X if (xx != a) {
X if (xx > a) {
X ++np;
X if (xx < xp) xp = xx;
X } else if (xx < a) {
X ++nm;
X if (xx > xm) xm = xx;
X }
X sum += dum=1.0/(eps+fabs((double)(xx-a)));
X sumx += xx * dum;
X }
X }
X }
X
X stemp = (sumx / sum) - a;
X if ((np-nm) >= 2) {
X /* guess is too low make another pass */
X am = a;
X aa = stemp < 0.0? xp : xp + stemp*AMP;
X if (aa > ap) aa = 0.5*(a+ap);
X eps = AFAC * fabs((double)(aa-a));
X a = aa;
X } else if ((nm-np) >= 2) {
X /* guess is too high make another pass */
X ap = a;
X aa = stemp > 0.0? xm : xm+stemp*AMP;
X if (aa < am) aa = 0.5*(a+am);
X eps = AFAC*fabs((double)(aa-a));
X a = aa;
X } else { /* got it */
X if (n%2 == 0) {
X Median = 0.5*(np==nm? xp+xm : np>nm? a+xp : xm+a);
X } else {
X Median = np == nm? a : np>nm? xp : xm;
X }
X return(Median);
X }
X }
X}
END_OF_FILE
if test 25452 -ne `wc -c <'Csrc/plotting.c'`; then
echo shar: \"'Csrc/plotting.c'\" unpacked with wrong size!
fi
chmod +x 'Csrc/plotting.c'
# end of 'Csrc/plotting.c'
fi
echo shar: End of archive 3 \(of 3\).
cp /dev/null ark3isdone
MISSING=""
for I in 1 2 3 ; do
if test ! -f ark${I}isdone ; then
MISSING="${MISSING} ${I}"
fi
done
if test "${MISSING}" = "" ; then
echo You have unpacked all 3 archives.
rm -f ark[1-9]isdone
else
echo You still need to unpack the following archives:
echo " " ${MISSING}
fi
## End of shell archive.
exit 0
--
Mail submissions (sources or binaries) to <amiga@uunet.uu.net>.
Mail comments to the moderator at <amiga-request@uunet.uu.net>.
Post requests for sources, and general discussion to comp.sys.amiga.