Back to index

radiance  4R0+20100331
gcalc.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: gcalc.c,v 1.3 2003/11/15 02:13:37 schorsch Exp $";
00003 #endif
00004 /*
00005  *  gcalc.c - routines to do calculations on graph files.
00006  *
00007  *     7/7/86
00008  *     Greg Ward Larson
00009  */
00010 
00011 #include  <stdio.h>
00012 #include  <string.h>
00013 
00014 #include  "mgvars.h"
00015 
00016 
00017 static double  xsum, xxsum, ysum, yysum, xysum;
00018 static double  xmin, xmax, ymin, ymax;
00019 static double  lastx, lasty, rsum;
00020 static int  npts;
00021 
00022 static void gcheader(char  *types);
00023 static void calcpoint(int  c, double  x, double  y);
00024 static void gcvalue(int  c, char  *types);
00025 
00026 void
00027 gcalc(               /* do a calculation */
00028        char  *types
00029 )
00030 {
00031        int  i;
00032 
00033        if (strchr(types, 'h') == NULL)
00034               gcheader(types);
00035        
00036        xmin = gparam[XMIN].flags & DEFINED ?
00037                      varvalue(gparam[XMIN].name) :
00038                      -1e10;
00039        xmax = gparam[XMAX].flags & DEFINED ?
00040                      varvalue(gparam[XMAX].name) :
00041                      1e10;
00042 
00043        for (i = 0; i < MAXCUR; i++) {
00044               xsum = xxsum = ysum = yysum = xysum = 0.0;
00045               rsum = 0.0;
00046               npts = 0;
00047               mgcurve(i, calcpoint);
00048               gcvalue(i, types);
00049        }
00050 }
00051 
00052 
00053 void
00054 gcheader(                   /* print header */
00055        register char  *types
00056 )
00057 {
00058        printf("__");
00059        while (*types)
00060               switch (*types++) {
00061               case 'n':
00062                      printf("|_Label___________");
00063                      break;
00064               case 'a':
00065                      printf("|____Mean______S.D._");
00066                      break;
00067               case 'm':
00068                      printf("|_____Min_______Max_");
00069                      break;
00070               case 'i':
00071                      printf("|___Integ_");
00072                      break;
00073               case 'l':
00074                      printf("|___Slope_____Intcp______Corr_");
00075                      break;
00076               default:
00077                      break;
00078               }
00079        printf("\n");
00080 }
00081 
00082 
00083 void
00084 gcvalue(             /* print the values for the curve */
00085        int  c,
00086        register char  *types
00087 )
00088 {
00089        double  d1, d2, d3, sqrt();
00090 
00091        if (npts < 1)
00092               return;
00093 
00094        printf("%c:", c+'A');
00095        while (*types)
00096               switch (*types++) {
00097               case 'n':
00098                      printf("  %-16s", cparam[c][CLABEL].flags & DEFINED ?
00099                                    cparam[c][CLABEL].v.s : "");
00100                      break;
00101               case 'a':
00102                      d1 = sqrt((yysum - ysum*ysum/npts)/(npts-1));
00103                      printf(" %9.2f %9.3f", ysum/npts, d1);
00104                      break;
00105               case 'm':
00106                      printf(" %9.2f %9.2f", ymin, ymax);
00107                      break;
00108               case 'i':
00109                      printf(" %9.2f", rsum);
00110                      break;
00111               case 'l':
00112                      d3 = xxsum - xsum*xsum/npts;
00113                      d1 = (xysum - xsum*ysum/npts)/d3;
00114                      d2 = (ysum - d1*xsum)/npts;
00115                      d3 = d1*sqrt(d3/(yysum - ysum*ysum/npts));
00116                      printf(" %9.5f %9.2f %9.5f", d1, d2, d3);
00117                      break;
00118               default:
00119                      break;
00120               }
00121        printf("\n");
00122 }
00123 
00124 
00125 void
00126 calcpoint(           /* add a point to our calculation */
00127        int  c,
00128        double  x,
00129        double  y
00130 )
00131 {
00132        if (x < xmin || x > xmax)
00133               return;
00134 
00135        xsum += x;
00136        xxsum += x*x;
00137        ysum += y;
00138        yysum += y*y;
00139        xysum += x*y;
00140 
00141        if (npts) {
00142               if (y < ymin)
00143                      ymin = y;
00144               else if (y > ymax)
00145                      ymax = y;
00146        } else
00147               ymin = ymax = y;
00148 
00149        if (npts)
00150               rsum += ( y + lasty )*( x - lastx )/2.0;
00151        lastx = x;
00152        lasty = y;
00153 
00154        npts++;
00155 }