root/trunk/EBaporate/ebaporate.c

/***************************************
  $Header$

  This file has the ebaporate front-end to the code, which inputs a PPM file
  with a hearth image and integrates a function over it, such as the
  evaporation rate or Al/Ti ratio.
***************************************/


#include <stdio.h>
#include <string.h>
#include <math.h>

#define COLS 1     /* X cell size */
#define ROWS 1     /* Y cell size */
#define THRESH 37  /* Minimum color value to be considered in Intens-T curve */
#define COLOR 2    /* Color being used by the function (1red, 2green, 3blue) */
#define POWER 1    /* Power of output, 2 accentuates high-temp differences */

/* Clausius-Clapeyron constants */
/* Al/Ti, for gamma calculation */
#define AA (16450-23200)
#define BB (12.36-11.74)
#define CC (-1.023+0.66)
#define DD 0.
#define CONST 1.33216 /* =sqrt(Mti/Mal) */
/* V/Ti, for gamma calculation
#define AA (26900-23200)
#define BB (10.12-11.74)
#define CC (0.33+0.66)
#define DD (-.265)
#define CONST 0.974292887 /* =sqrt (Mv/Mal) */

#define PI 3.141592653589793 /* 'Nuff said */


/*++++++++++++++++++++++++++++++++++++++
  The usual main() function.

  int argc Argument count.

  char *argv[] Argument values.
  ++++++++++++++++++++++++++++++++++++++*/

void main (int argc, char *argv[])
{
  int i, j, n=0, rowd, cold;
  double *rower, total=0., last, func(), fTH, f255;
  unsigned char insup [80], *outs;
  FILE *ppmout;

  /*+ Get the top of the file, verify PPM, and determine rows and columns. +*/
  scanf ("%s", insup);
  if (strcmp (insup, "P6"))
    exit (printf ("Not PPM.\n"), 1);

  scanf ("%s", insup);
  if (insup [0] == '#')
    {
      while (insup [1] != 13 && insup [1] != 10)
        scanf ("%c", insup+1);
      scanf ("%s", insup);
    }

  sscanf (insup, "%d", &cold);
  scanf ("%d", &rowd);
  scanf ("%d", &i);
  printf ("Picture rows: %d, columns: %d, bytemask:%d\n", rowd, cold, i);
  scanf ("%c", insup+1);

  /*+ Open the output temperature map file. +*/
  if ((ppmout = fopen ("map.ppm", "w")) == NULL)
    exit (printf ("No ppm file!\n"), 1);
  fprintf (ppmout, "P5\n# CREATOR: Adam Powell's ebaporate software\n");
  fprintf (ppmout, "%d %d\n255\n",
          (cold+COLS-1)/COLS, (rowd+ROWS-1)/ROWS);

  /*+ Allocate memory for everything. +*/
  if ((rower = (double *)malloc (sizeof (double)*((cold+COLS-1)/COLS))) == NULL)
    exit (printf ("Not enough memory...\n"));
  if ((outs = (char *)malloc (sizeof (char)*((cold+COLS-1)/COLS))) == NULL)
    exit (printf ("Not enough memory...\n"));
  for (j=0; j<(cold+COLS-1)/COLS; j++)
    rower [j]=0.;

  /*+ The big loop: read the input PPM file, calculate the local function value
     and add it to the total, and store the output PPM. +*/
  for (i=0; i<rowd; i++)
    {
      for (j=0; j<cold; j++)
        {
          scanf ("%c%c%c", insup+1, insup+2, insup+3);
          if (insup [COLOR] < THRESH)
            rower [j/COLS] = 256. * ROWS * COLS;
          rower [j/COLS] += insup [COLOR];
        }
      fTH = func ((double)THRESH);
      f255=func (255.);
      if (i%ROWS == ROWS-1 || i == rowd-1)
        for (j=0;j<(cold+COLS-1)/COLS;j++)
          {
            last = rower [j] /= ROWS*COLS;
            if (rower [j] < 256)
              {
                total += func (rower [j]);
                n++;
                outs [j] = (unsigned char)
                  (220. - 220.*pow ((func (rower [j])-fTH)/(f255-fTH), POWER));
              }
            else
              {
                outs [j] = 255;
              }
            fprintf (ppmout, "%c", outs [j]);
            rower [j] = 0.;
          }
    }
  fclose (ppmout);

  /* Make temperature bar ppm */
  printf ("Making Tbar...\n");
  if ((ppmout = fopen ("tbar.ppm", "w")) == NULL)
    exit (printf ("No Tbar file!\n"), 1);
  fprintf (ppmout, "P5\n# CREATOR: Adam's evapometer software\n216 20\n255\n");
  for (i=0; i<20;i++)
    {
      for (j=0; j<24;j++)
        fprintf (ppmout, "%c", 255);
      for (j=24; j<216;j++)
        fprintf (ppmout, "%c", (unsigned char)
                 (220.-220.*pow(((double)j-24.)/192., POWER)));
    }
  fclose (ppmout);

  printf ("Min, max function values: %lf %lf\n", fTH, f255);
  printf ("(Last rower: %.2lf; last chars: %d %d %d)\n", last,
          *(insup+1), *(insup+2), *(insup+3));
  printf ("Average functional value: %le, based on %d out of %d cells.\n\n",
          total/n, n, ((cold+COLS-1)/COLS)*((rowd+ROWS-1)/ROWS));
}


/*++++++++++++++++++++++++++++++++++++++
  This function converts color values (0-255) into aluminum activity
  coefficients.  It can easily be modified to give evaporation rates, and at
  some point there will be a function pointer to allow multiple functions to
  handle this in different ways.

  double func It returns evaporation rate, activity coefficient, etc.

  double color Color value to evaluate.
  ++++++++++++++++++++++++++++++++++++++*/

double func (double color)
{
  double temp, vp;
  temp = 288.31 * log (color*23.02927);          /* Temperature, green-based */
  /* temp = 322.58 * log (color*3.0794)+273;       /* Temperature, red-based */
  vp = exp (log (10.) * (-AA/temp + BB + CC*log10 (temp) + DD*.001*temp));
  /* Vap press */
  vp *= CONST;
  /*vp *= .04788*101300./760./sqrt (2*PI*0.04788*8.314*temp); /* 101300./760.*6.2213647e-07/sqrt (temp); */
  return (temp); }