root/trunk/EBaporate/ebsurf.c

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

  This is the brains of the dynamic evaporation estimator.  It is linked into
  the ebsurfweb and ebaporate binaries.

  It does a multi-cycle 1-D transient simulation of heat conduction down into
  the melt as the beam scans over a spot on the surface, as described in:

  A. Powell, J. Van Den Avyle, B. Damkroger, J. Szekely and U. Pal
  ``Analysis of Multicomponent Evaporation in Electron Beam Melting and
  Refining of Titanium Alloys,''
  +latex+{\em Metall. Trans.} {\bf 38B}, 1227-1239 (1997).
  -latex-Metall. Trans. 38B, 1227-1239 (1997).

  Features include:
  +latex+\begin{itemize} \item
  +html+ <ul><li>
  Multiple postscript graph outputs.
  +latex+\item
  +html+ </li><li>
  Beam can hit once per cycle (default) or follow programmable patterns.
  +latex+\item
  +html+ </li><li>
  Variable timestep size to accurately capture beam strike and efficiently
  calculate relaxation.
  +latex+\item
  +html+ </li><li>
  Calculation of beam penetration depth from voltage, using function given
  by Schiller (see shorts.c for details).
  +latex+\item
  +html+ </li><li>
  Volumetric heat generation due to beam through the penetration depth,
  using functionality given by Schiller.
  +latex+\item
  +html+ </li><li>
  Crank-Nicholson timestepping with Newton iteration for nonlinear boundary
  conditions.
  +latex+\end{itemize}
  +html+ </ul>

  This program uses cgs units and Kelvin.
***************************************/


#include "ebsurf.h"


/*++++++++++++++++++++++++++++++++++++++
  This function runs through multiple time-periodic cycles of timesteps until
  steady-state is reached for the base temperature or surface temperature T0
  (depending on whether the TSUR flag is set in m).  ``Steady-state'' is
  determined by the difference between surface temperature at the start and end
  of the cycle.

  int cycles It returns zero or an error code.

  double *X Pointer to storage for positions (depth below surface).

  double *t Pointer to storage for time step times.

  double *T Pointer to storage for temperature array.

  double T0 Base, or average surface, temperature.

  double *big Pointer to storage for tridiagonal matrix coefficients.

  struct param *p EBsurf Parameter structure.

  int m Flags: TSUR if T0 is to be the surface temperature, OUTPUT to print
  graphs, GREET to print other information such as timestep completion.

  double *ou Little array of output values.

  FILE *fule File name for HTML template.

  int pt Increment to use in ou array.
  ++++++++++++++++++++++++++++++++++++++*/

int cycles (double *X, double *t, double *T, double T0, double *big,
            struct param *p, int m, double *ou, FILE *fule, int pt)
{
  double Gto, tmp;
  double timestep(), radloss(), loss(), vloss(), vloss2(), gtnorm(), gtdisk();
  int i=0,j=0,s;
  char *a,inp[80];
  void intcond(), output(), printto(), skipto();

  /*+ First we call initcond() to fill in the depth (X), time (t), and initial
    temperature (T) arrays, along with the generation (G) and tridiagonal
    matrix arrays (A,B,C) stored in ``big'', etc. +*/
  initcond (X, t, T, T0, (m&TSUR)?0.:p->Tsurf, big, p);
  if (m & TSUR)
    T0=T[0];

  /*+ The big loop then runs the timesteps.  This can happen in one of two
    ways: +*/
  Gto = p->gtime (0., p);
  if (p->cyc)
    {
      /*+ If the "cyc" member of the param structure p is nonzero, then it will
        run through cyc beam scan cycles, storing all of those cycles for
        post-processing. +*/
      for (j=0;j<p->cyc;j++)
        {
          printf ("Cycle %d",j+1);
          if (m)
            printf (", timesteps:");
          for (i=0;i<p->ts;i++)
            Gto=timestep(*p,X,t,T+(p->zs+1)*i,big,Gto,i+1,m);
          printf (" dT=%.3lf, Tnew=%.3lf\n", T[(p->zs+1)*p->ts] - T[0],
                  T[(p->zs+1)*p->ts]);
          T += p->ts*(p->zs+1);
        }
      T -= p->cyc*p->ts*(p->zs+1);
    }
  else
    {
      /*+ Otherwise, it will run through as many beam scan cycles as it takes
        to converge the surface temperature to steady-state, storing only the
        last cycle for post-processing. +*/
      j=0;
      do
        {
          if(j)
            {
              /* Calculate the "surface average temperature". */
              tmp = (T[p->zs-1] * t[1] +
                     T[(p->zs+1)*(p->ts+1)-2] * (p->tfin - t[p->ts-1])) * 0.5;
              for (i=2; i<=p->ts; i++)
                tmp += T[(p->zs+1)*i-2] * (t[i] - t[i-2])*0.5;
              tmp = (tmp/p->tfin-T[p->zs])*p->w/(p->w-X[p->zs-1]) + T[p->zs];
              if (m&TSUR)
                {     /* Move 1st-step Ts so it's right */
                  for (i=0; i<=p->zs; i++)
                    T[i] = T[p->ts*(p->zs+1)+i] + T0-tmp;
                }
              else
                {
                  for (i=0; i<=p->zs; i++)
                    T[i] = T[p->ts*(p->zs+1)+i];
                }
              printf ("<br>After cycle %d, surface average temperature=%g<sup>o</sup>C\n",
                     j, tmp-273.15);
            }
          for (i=0; i<p->ts; i++)
            {
              Gto = timestep (*p, X, t, T+(p->zs+1)*i, big, Gto, i+1, m);
            }
        }
      while(++j<MAXCYC && _ABS(T[0] - T[(p->zs+1)*p->ts]) > p->ct);
    }

  if (!p->cyc && _ABS(T[0] - T[(p->zs+1)*p->ts]) <= p->ct)
    printf("<p><b>Converged</b> after ");
  else
    printf("Completed ");
  printf("%d cycles.<p>",j);

  s=(p->cyc?p->cyc-1:0)*p->ts;
  if(m&TSUR)
    T0= *(T+p->zs);
  if(p->cyc)
    T-=j*p->ts*(p->zs+1);
  i=p->cyc;
  if(p->cyc==0)
    p->cyc=1;

#define DELTAT ((j?*(t+j+1)-*(t+j-1):*(t+1)+p->tfin-*(t+p->ts-1))*0.5)
  /* Calculate average temperature, evap (avg temp), evaporation */
  ou[0] = ou[2*pt] = 0.;
  for(j=0;j<p->ts;j++)
    {
      ou[0] += T[(s+j)*(p->zs+1)+p->zs-1] * DELTAT;
      /* Replace vloss here with vloss2 for aluminum k'' calculation */
      ou[2*pt] += vloss (T[(s+j)*(p->zs+1)], p->E, p->A, p->C, p->D) * DELTAT;
    }
  ou[0] = (ou[0]/p->tfin - T0) *p->w/(p->w-X[p->zs-1])+T0; /* p->Ma*10.* */
  ou[2*pt] /= p->tfin;
  ou[pt] = vloss(ou[0], p->E, p->A, p->C, p->D); /* here too! */

  /* Report: results table: Tsurf, const T ev rate, real ev rate... */
  printto (fule,"<TD>1</TD>");
  printf ("<TD><A HREF=\"http://lyre.mit.edu/~powell/Software/outputs.html#tsurf\"><I>T</I><SUB>surf</SUB></A>=%g<SUP>o</SUP>C\n",ou[0] - 273.15);
  printto (fule,"<TD>1</TD>"); printf("<TD>%g</TD>\n", ou[pt] * p->Ma*1e4);
  printto (fule,"<TD>1</TD>"); printf("<TD>%g</TD>\n", ou[2*pt] * p->Ma*1e4);
  printto (fule,"<TD>1</TD>");
  printf ("<TD>%g%%</TD>\n", (ou[2*pt] - ou[pt])*100. / ou[pt]);

  /* ..CT heat loss, real heat loss, CT & real rad loss, CT & real ev loss.. */
  printto (fule,"<TD>1</TD>"); printf("<TD>%g</TD>\n",loss(*ou,p)/1.e7);
  tmp=0.;
  for (j=0; j<p->ts; j++)
    tmp += loss (T[(s+j)*(p->zs+1)], p) * DELTAT;
  printto (fule,"<TD>1</TD>"); printf("<TD>%g</TD>\n",tmp/p->tfin/1.e7);
  printto (fule,"<TD>1</TD>");
  printf ("<TD>%g%%</TD>\n",(tmp/p->tfin-loss(*ou,p))*100./loss(*ou,p));

  printto (fule,"<TD>1</TD>"); printf("<TD>%g</TD>\n",radloss(*ou,p)/1.e7);
  tmp=0.; for(j=0;j<p->ts;j++) tmp+=radloss(*(T+(s+j)*(p->zs+1)),p)*DELTAT;
  printto(fule,"<TD>1</TD>"); printf("<TD>%g</TD>\n",tmp/p->tfin/1.e7);
  printto(fule,"<TD>1</TD>");
  printf("<TD>%g%%</TD>\n",(tmp/p->tfin-radloss(*ou,p))*100./radloss(*ou,p));

  printto(fule,"<TD>1</TD>"); printf("<TD>%g</TD>\n",p->Q0*exp(p->S* *ou)*
                                     vloss(*ou,p->E,p->A,p->C,p->D)/1.e7);
  tmp=0.; for(j=0;j<p->ts;j++)
    tmp+=p->Q0*exp(p->S* *(T+(s+j)*(p->zs+1)))*
      vloss(*(T+(s+j)*(p->zs+1)),p->E,p->A,p->C,p->D)*DELTAT;
  printto(fule,"<TD>1</TD>"); printf("<TD>%g</TD>\n",tmp/p->tfin/1.e7);
  printto(fule,"<TD>1</TD>");
  printf("<TD>%g%%</TD>\n",(tmp/p->tfin-
                          p->Q0*exp(p->S* *ou)*vloss(*ou,p->E,p->A,p->C,p->D))
         *100./(p->Q0*exp(p->S* *ou)*vloss(*ou,p->E,p->A,p->C,p->D)));

  /* ...beam heat input, layer bottom temperature, conduction at bot, surf. */
  printto(fule,"<TD>1</TD>"); /* Theoretical centerline heat flux */
  printf("<TD>%g</TD>\n",p->pow*4./M_PI/p->spx/p->spy/1.e10);
  printto(fule,"<TD>1</TD>"); /* Actual centerline heat flux = integral of */
  tmp=*big; for(j=1;j<p->zs;j++) tmp+= *(big+j); /* q(z) dz */
  printf("<TD>%g</TD>\n",tmp*p->k/p->al/1.e10);
  printto(fule,"<TD>1</TD>"); if(p->gtime==gtnorm)
    printf("<TD>%g</TD>\n", /* Theoretical beam heat input (sqrt(PI)/2) */
           p->pow*2./sqrt(M_PI)/p->spx/p->spy*p->tres/p->tfin/1e7);
  else printf("<TD>%g</TD>\n", /* Even distribution for liss & disk */
              p->pow/p->plen/p->pwid/1e7*((p->gtime==gtdisk)?4./M_PI:1.));
  tmp=0.; for(j=0;j<p->ts;j++) tmp+=(p->gtime)(*(t+j),p)*DELTAT;
  printto(fule,"<TD>1</TD>"); /* Actual beam heat input */
  printf("<TD>%g</TD>\n",p->pow*4/M_PI/p->spx/p->spy*tmp/p->tfin/1.e7);
  printto(fule,"<TD>1</TD>"); printf("<TD>%g</TD>\n",T0-273.15); /* Layer bot temp */
  tmp=0.; for(j=0;j<p->ts;j++) tmp+= *(T+(s+j+1)*(p->zs+1)-2)*DELTAT;
  printto(fule,"<TD>1</TD>"); /* Conduction through bottom */
  printf("<TD>%g</TD>\n",(tmp/p->tfin-T0)*p->k/(p->w-*(X+p->zs-1))/1.e7);
  tmp=0.; for(j=0;j<p->ts;j++)
    tmp+= (*(T+(s+j)*(p->zs+1))-*(T+(s+j)*(p->zs+1)+1))*DELTAT;
  printto(fule,"<TD>1</TD>"); /* Approx. conduction through surface */
  printf("<TD>%g</TD>\n",p->k*tmp/ *(X+1)/p->tfin/1.e7);
#undef DELTAT

  /* Output graphs (assumes last cycle is always output) */
  if((m&OUTPUT)!=0) {
    p->Tin=T0; printto(fule,"<TD>1</TD>");
    if(!p->pout) skipto(fule,"<TD>1</TD>"); else {
      output(p,X,t,big,PATT); printto(fule,"<TD>1</TD>"); } /* Pattern */
    if(!p->gout) skipto(fule,"<TD>1</TD>"); else {
      output(p,X,t,big,GENZ); printto(fule,"<TD>1</TD>"); } /* Heat gen func */
    if(!p->gtou) skipto(fule,"<TD>1</TD>"); else {
      output(p,X,t,big,GENT); printto(fule,"<TD>1</TD>"); } /* Heat gen hist */
    if(!p->fout) skipto(fule,"<TD>1</TD>"); else {
      output(p,X,t,T,FLUX); printto(fule,"<TD>1</TD>"); }   /* Bottom fluxes */
    if(!p->sout) skipto(fule,"<TD>1</TD>"); else {
      output(p,X,t,T,SURF); printto(fule,"<TD>1</TD>"); }  /* All surf temps */
    if(!p->tout) skipto(fule,"<TD>1</TD>"); else {
      output(p,X,t,T,TMPS); printto(fule,"<TD>1</TD>"); } /* All temps (3-D) */
    output(p,X,t,T+(i?(i-1)*(p->zs+1)*(p->ts):0),ONEC);    /* Last cyc temps */
  }
  p->cyc=i;

  return 0;
}


/*++++++++++++++++++++++++++++++++++++++
  Perform one timestep in a cycle.

  double timestep It returns a double equal to the heat generation per unit
  area at the next timestep.

  struct param p EBsurf parameter structure.

  double *X The array of z-values.

  double *t The array of time step times.

  double *T The big array of temperatures, but pointing to this timestep.

  double *big Array storing lots of different z-dependent parameters.

  double Gto Heat generation rate per unit area at this timestep.

  int s Timestep number for indexing into t array.

  int m Flags, in particular GREET is important here.
  ++++++++++++++++++++++++++++++++++++++*/

double timestep(struct param p, double *X, double *t, double *T, double *big,
                double Gto, int s, int m)
{
  int i;
  void tridag();
  double Gt,dt,Tb,*Gx,*A,*Ag,*At,*B,*Bg,*Bt,*C,*Cg,*Ct,*O,*R,loss(),lossder();

  Gx = big;
  A  = big + p.zs;
  Ag = big + p.zs*2;
  At = big + p.zs*3;
  B  = big + p.zs*4;
  Bg = big + p.zs*5;
  Bt = big + p.zs*6;
  C  = big + p.zs*7;
  Cg = big + p.zs*8;
  Ct = big + p.zs*9;
  O  = big + p.zs*10;
  R  = big + p.zs*11;

  /*+ First it initializes the parameters, particularly the tridiagonal matrix,
    according to the size of this timestep. +*/
  dt = 1./(t[s]-t[s-1]);
  Tb = (1.-p.th)/p.th;
  Gt = (p.gtime)(t[s], &p);
  for(i=0;i<p.zs;i++)
    {
      A[i] = Ag[i] + At[i]*dt;
      B[i] = Bg[i] + Bt[i]*dt;
      C[i] = Cg[i] + Ct[i]*dt;
      O[i] = (Ag[i]*Tb - At[i]*dt) * T[i-1] + (Bg[i]*Tb - Bt[i]*dt)* T[i]+
        (Cg[i]*Tb - Ct[i]*dt) * T[i+1] - ((1-p.th)*Gto+p.th*Gt) * Gx[i];
      T[p.zs+1+i] = T[i];
    }

  if(m&GREET)
    printf(" %d",s);

  O[0] += p.al * (1-p.th) * loss (T[0], &p)/p.k;
  T[2*p.zs+1] = T[p.zs];
  T+=p.zs+1;
  s=0;

  /*+ Then it uses Newton-Raphson iteration to perform one semi-implicit
    timestep, finishing when the difference between first and last temperature
    is less than the "it" member of the param structure p. +*/
  do
    {
      Tb = T[0];
      B[0] = Bg[0] + Bt[0]*dt;
      for(i=0;i<p.zs;i++)
        R[i] = O[i] + (i ? (A[i]*T[i-1]) : 0) + B[i]*T[i] + C[i]*T[i+1];
      R[0] += p.al * p.th * loss (T[0], &p)/p.k;
      B[0] += p.al * p.th * lossder(T[0], &p)/p.k;
      tridag (A, B, C, R, T, p.zs);
    }
  while(_ABS(*T-Tb)>p.it);

  return Gt;
}


/*++++++++++++++++++++++++++++++++++++++
  Tridiagonal matrix solver, solves (a,b,c)arr=r, and u+=arr.

  double *a Subdiagonal in tridiagonal matrix.

  double *b Diagonal in tridiagonal matrix.

  double *c Superdiagonal in tridiagonal matrix.

  double *r Right-hand side.

  double *u Array to add the solution to.

  int n Number of members in the array.
  ++++++++++++++++++++++++++++++++++++++*/

void tridag(double *a, double *b, double *c, double *r, double *u, int n)
{
  int j;
  double bet,*gam,*arr;

  if (!(gam=MALLOC(double,2*n)))
    exit (printf ("No memory for gamma vector.\n"));
  arr=gam+n;

  if (b[0] == 0.0)
    exit (printf ("Matrix has a zero corner.\n"));

  arr[0] = r[0]/(bet=b[0]);
  for (j=1; j<n; j++)
    {
      gam[j] = c[j-1]/bet;
      bet = b[j] - a[j]*gam[j];
      if (bet == 0.0)
        exit( printf("Matrix has a zero pivot.\n"));
      arr[j] = (r[j] - a[j]*arr[j-1]) / bet;
    }

  u[n-1] -= arr[n-1];
  for (j=(n-2); j>=0; j--)
    u[j] -= arr[j] -= gam[j+1]*arr[j+1];

  free(gam);
}


/*++++++++++++++++++++++++++++++++++++++
  Calculate the initial conditions, generation function, tridiagonal matrix
  coefficients, etc.

  double *X The array of z-values.

  double *t The array of time step times.

  double *T The big array of temperatures.

  double T0 Initial base temperature.

  double TS Initial surface temperature.

  double *big Array storing lots of different z-dependent parameters.

  struct param *p EBsurf parameter structure.
  ++++++++++++++++++++++++++++++++++++++*/

void initcond (double *X, double *t, double *T, double T0, double TS,
               double *big, struct param *p)
{
  int i=0;
  double tmp,*Gx,*Ag,*At,*Bg,*Bt,*Cg,*Ct;

  /* Position of vectors within "big" */
  Gx=big;
  Ag=big+p->zs*2;
  At=big+p->zs*3;
  Bg=big+p->zs*5;
  Bt=big+p->zs*6;
  Cg=big+p->zs*8;
  Ct=big+p->zs*9;

  /* Make space and time grids */
  X[0] = 0.;
  X[p->zs] = p->w;
  t[0] = 0.;
  tmp = log(p->zr)/(p->zs-1);

  for (i=1; i<=p->zs; i++)
    X[i] = X[i-1] + exp(tmp*i);
  tmp = p->w / X[p->zs];
  for (i=1; i<p->zs; i++)
    X[i] *= tmp;
  X[p->zs] = p->w;

  tmp = log (p->tr) / (p->ts-1);
  for (i=1; i<=p->ts; i++)
    t[i] = t[i-1] + exp(tmp*i);
  tmp = p->tfin/t[p->ts];
  for(i=1; i<p->ts; i++)
    t[i] *= tmp;
  t[p->ts] = p->tfin;

  printf("Space  grid: ");
  for(i=0;i<=p->zs;i++)
    printf("%.5lf ", X[i]);
  printf("\nTime  grid: ");
  for(i=0;i<=p->ts;i++)
    printf("%.5lf ", t[i]);
  printf("\n");

  /* tmp = flux/k = -dT/dx; make initial condition with this slope */
  tmp = p->pow*4./p->k/M_PI/p->spx/p->spy; /* flux/k */
  if (TS!=0.)
    for (i=0;i<=p->zs;i++)
      T[i] = TS + X[i] / p->w * (T0-TS);
  else
    {
      for(i=0;i<=p->zs;i++)
        T[i] = T0 + tmp * p->tres/p->tfin * sqrt (M_PI)/2. * (p->w - X[i]);
    }

  /* Calculate matrices: Ag,Bg,Cg, At,Bt,Ct (alpha*theta*stiffness, mass) */
  for(i=1;i<p->zs;i++)
    {
      Ag[i] = -p->al*p->th/(*(X+i)-*(X+i-1));
      At[i] = (X[i] - X[i-1]) / 6;
      Bt[i] = (X[i+1] - X[i-1]) / 3;
      Bg[i] = Bt[i] * 3 * Ag[i] / (X[i] - X[i+1]);
      Cg[i] = -p->al*p->th/(X[i+1] - X[i]);
      Ct[i] = (X[i+1] - X[i])/6;
    }
  Ag[0] = At[0] = 0;
  Bg[0] = p->al*p->th / X[1];
  Bt[0] = X[1] / 3;
  Cg[0] = -Bg[0];
  Ct[0] = Bt[0] / 2;

  /* Calculate z-dependent generation function (tmp=flux/k from above).
     This calculates the integral of the actual Gaussian-fit function q(z) over
     the shapefunctions in each element phi_i(z), and scales that element's
     generation function accordingly. */
  printf ("Penetration depth=%g cm; Z1=%g cm; ratio=%g\n", p->dp, X[1],
          p->dp / X[1]);
  if (p->dp==0.)
    {
      *Gx=tmp*p->al;
      for(i=1;i<p->zs;i++)
        *(Gx+i)=0.;
    }
  else
    {
      Gx[0] = 0.;
      tmp *= 1.364242002*p->al/p->dp;
      for (i=0; i<p->zs; i++)
        {
#define Z0 (*(X+i-1))
#define Z1 (*(X+i))
          if(i)
            Gx[i]= -tmp/(Z1-Z0)*p->dp*p->dp/8.* /*Integral over prev interval*/
              (exp(-4*(Z1/p->dp-1./3.) * (Z1/p->dp-1./3.))-/*of phi_i(z)q(z)dz*/
               exp(-4*(Z0/p->dp-1./3.) * (Z0/p->dp-1./3.)))+
              tmp*p->dp/2.*(-Z0/(Z1-Z0)+p->dp/(Z1-Z0)/3.)* /*Erf error correct*/
              sqrt(M_PI)/2.*(erf(2.*Z1/p->dp-2./3.)-erf(2.*Z0/p->dp-2./3.));
#undef Z1
#undef Z0
#define Z0 (*(X+i))
#define Z1 (*(X+i+1))
          Gx[i] += tmp/(Z1-Z0)*p->dp*p->dp/8.*  /*Integral over next interval*/
            (exp(-4*(Z1/p->dp-1./3.) * (Z1/p->dp-1./3.))- /*of phi_i(z)q(z)dz*/
             exp(-4*(Z0/p->dp-1./3.) * (Z0/p->dp-1./3.)))+
            tmp*p->dp/2.*(Z1/(Z1-Z0)-p->dp/(Z1-Z0)/3.)*   /*Erf error correct*/
            sqrt(M_PI)/2.*(erf(2.*Z1/p->dp-2./3.)-erf(2.*Z0/p->dp-2./3.));
#undef Z1
#undef Z0
        }
    }
}


/*++++++++++++++++++++++++++++++++++++++
  Little replacement for strcmp (forgot what was wrong with the old one).

  int mystrcmp Returns 0 if equal, something else otherwise.

  char *string1 First string to compare.

  char *string2 Second string to compare.
  ++++++++++++++++++++++++++++++++++++++*/

int mystrcmp (char *string1, char *string2)
{
  int i=0;
  while (i<MIN(strlen(string1),strlen(string2)) && string2[i]==string1[i])
    i++;
  return string2[i] && string1[i];
}


/*++++++++++++++++++++++++++++++++++++++
  Ancillary function for paraminp, scans ``field'' string for ``search'' key
  string presence and returns the int following it.

  int inputint Returns the integer inputted.

  char *field String to search.

  char *search Key to search for.

  int ret Default value to return if nothing is found.
  ++++++++++++++++++++++++++++++++++++++*/

int inputint (char *field, char *search, int ret)
{
  int i;
  if (field)
    for (i=0; i+strlen(search) <= strlen(field); i++)
      if (!(mystrcmp (field+i, search)))
        sscanf (field+strlen(search)+i, "%d", &ret);
  return ret;
}


/*++++++++++++++++++++++++++++++++++++++
  Ancillary function for paraminp, scans ``field'' string for ``search'' key
  string presence and returns the floating point number following it.

  double inputdouble Returns the floating point number in double form.

  char *field String to search.

  char *search Key to search for.

  double ret Default value to return if nothing is found.
  ++++++++++++++++++++++++++++++++++++++*/

double inputdouble (char *field, char *search, double ret)
{
  int i;
  if (field)
    for (i=0; i+strlen(search) <= strlen(field); i++)
      if (!(mystrcmp (field+i, search)))
        sscanf (field+strlen(search)+i, "%lf", &ret);
  return ret;
}


/*++++++++++++++++++++++++++++++++++++++
  Scan qstring for parameters and store them in ebsurf parameter structure p,
  using defaults for those not found.

  char *paraminp It returns NULL.

  struct param *p Parameter structure in which to place parameters.

  char *qstring String to scan for parameters.
  ++++++++++++++++++++++++++++++++++++++*/

char *paraminp(struct param *p, char *qstring)
{
  int i=0,j=0;
  double a,t,t2;
  FILE *status;
  static char gfile[80],hfile[80],tfile[80],sfile[80],lfile[80],ffile[80],
  pfile[80];

  /* Material properties, defaults for liquid titanium */
  p->al=inputdouble(qstring,"alpha=", 0.1);
  p->k=inputdouble(qstring,"cond=", 30.)*1e5;
  p->Ma=inputdouble(qstring,"molmas=", 47.88);
  p->eps=inputdouble(qstring,"epsilon=", 0.5);            /* Rad. emissivity */
  p->Tenv=inputdouble(qstring,"envtemp=",20.)+273.15;    /* Cels->internal K */
  p->rho=inputdouble(qstring,"density=", 4.1);
  p->mu=inputdouble(qstring,"viscos=", 0.052);
  p->dsig=inputdouble(qstring,"dsigdT=", 2.6e-4);
  p->Q0=inputdouble(qstring,"q0vap=", 482.)*1e10;   /* Q0 in kJ/mol->erg/mol */
  p->S=-inputdouble(qstring,"qSvap=", 4.64e-5);       /* S in inverse Kelvin */
  p->A=inputdouble(qstring,"ccA=",23200.)*log(10.);    /* Claus-Clap log->ln */
  p->B=inputdouble(qstring,"ccB=", 11.74)*log(10.);
  p->C=inputdouble(qstring,"ccC=",-0.66)-0.5;   /* -0.5 for the sqrt in lang */
  p->D=inputdouble(qstring,"ccD=",0.)*log(10.)/1000.;
  p->E=TORR2CGSP*exp(p->B)/sqrt(2*M_PI*GR*p->Ma);          /* Langmuir coeff */

  /* Operating parameters */
  p->pow=inputdouble(qstring,"power=", 200.)*1e10;   /* Net power, kW->erg/s */
  p->volt=inputdouble(qstring,"voltage=", 30.);                        /* kV */
  if((p->dp=inputdouble(qstring,"pendep=", 0.)/10000.)==0.)
    p->dp=depvolt(p->volt)/p->rho; /* Microns penetr, default calc from volt */
  p->tfin=inputdouble(qstring,"mfreq=", 200.); /* for now tfin has frequency */
  p->xfrq=inputdouble(qstring,"xfreq=", 6.);   /* X and y freqs for patterns */
  p->yfrq=inputdouble(qstring,"yfreq=", 7.);
  if(p->tfin==0.)
    errer("Zero frequency not permitted.");
  p->plen=inputdouble(qstring,"patlen=", 240.);       /* Beam pattern length */
  p->pwid=inputdouble(qstring,"patwid=", 7.5);         /* Width for patterns */
  p->x=inputdouble(qstring,"locX=", 0.);       /* X and y where we calculate */
  p->y=inputdouble(qstring,"locY=", 0.);
  p->vb=inputdouble(qstring,"vbeam=",   /* Default velocity=patlen*frequency */
                    p->plen*p->tfin);
  if(p->vb==0.)
    errer("Zero beam velocity and pattern length not permitted.");
  p->tfin=1./p->tfin;                               /* Now tfin holds period */
  if((p->spx=inputdouble(qstring,"spsize=", 2.5))==0.||
     (p->spy=inputdouble(qstring,"spwid=", 2.5))==0.)
    errer("Zero spot area not permitted.");
  p->tres=inputdouble(qstring,"dwell=",     /* Default restime=size/velocity */
                         p->spx/p->vb);
  p->Tsurf=inputdouble(qstring,"Tsurf=", 1750.);
  p->Tin=inputdouble(qstring,"Tbot=", 0.);
  p->Tsurf += (p->Tsurf==0.) ? 0. : 273.15;           /* Celsius->internal K */
  p->Tin += (p->Tin==0.) ? 0. : 273.15;

  /* Simulation parameters */
  p->zs=inputint(qstring,"zsteps=", 45);             /* Steps in z-direction */
  p->zr=inputdouble(qstring,"zrat=", 1.);            /* Ratio max/min deltaz */
  p->ap=inputdouble(qstring,"aleph=", 1.);/* Used for default timestep ratio */
  p->bp=inputdouble(qstring,"bapheth=", 9.);   /* Used for default thickness */
  p->w=inputdouble(qstring,"width=",         /* Thickness of layer simulated */
                   p->dp + sqrt(p->bp*p->al*p->tfin));
  p->b=inputdouble(qstring,"bsplit=", 5.);  /* Approx. timesteps during tres */
  if((p->tr=inputdouble(qstring,"trat=",0.))==0.)    /* Ratio max/min deltat */
    {
      if (p->zr==1) /* calculate t, the smallest element size */
        t=p->w/p->zs; /* straightforward if z ratio = 1 */
      else
        {
          a=exp(log(p->zr)/(p->zs-1)); /* ratio of adjacent element sizes */
          if(a<1.)
            a=1./a; t=0;
          for(j=0;j<p->zs;j++)
            t+=pow(a,(double)j);
          t=p->w/t;
        }
      if(p->tres/p->b > t*t/p->al/p->ap)
        {
          p->tr=1; /* stability criterion/aleph smaller than bsplit/tres */
          if(p->ts==0)
            p->ts=p->tfin/t/t*p->al*p->ap;
        }
      else
        p->tr= (t*t/p->al/p->ap) / (p->tres/p->b);
    }
  if ((p->ts=inputint(qstring,"tsteps=", 0))==0)
    {
      t=1.;
      for (p->ts=2; p->tfin/t>p->tres/p->b; p->ts*=2)
        {
          a=exp(log(p->tr)/(p->ts-1));
          if (a<1.)
            a=1./a;
          t=0;
          for (j=0;j<p->ts;j++)
            t+=pow(a,(double)j);
        }
      for (i=p->ts/2; i>1; i/=2)
        {
          if(p->tfin/t>p->tres/p->b)
            p->ts+=i;
          else
            p->ts-=i;
          t=0;
          a=exp(log(p->tr)/(p->ts-1));
          for(j=0;j<p->ts;j++)
            t+=pow(a,(double)j);
        }
      if(p->tfin/t>p->tres/p->b) p->ts++;
    }
  p->th=inputdouble(qstring,"theta=",0.5);   /* Implicitness, 0.5=Crank-Nich */
  p->it=inputdouble(qstring,"itol=",0.0001); /* Max dT for timestep converge */
  p->ct=inputdouble(qstring,"ctol=",0.01);     /* Max dT for cycles converge */
  p->cyc=inputint(qstring,"cyc=", 0);/* Number of cycles to run, 0->converge */

  /* Calculate real Aleph and Bapheth, in case tr and ts were set. */
  if(p->zr==1)
    t=p->w/p->zs;
  else
    {
      a=exp(log(p->zr)/(p->zs-1));
      if(a<1.)
        a=1./a;
      t=0;
      for(j=0;j<p->zs;j++)
        t+=pow(a,(double)j);
      t=p->w/t;
    }
  if(p->tr==1)
    t2=p->tfin/p->ts;
  else
    {
      a=exp(log(p->tr)/(p->ts-1));
      if(a>0.)
        a=1./a;
      t2=0;
      for(j=0;j<p->ts;j++)
        t2+=pow(a,(double)j);
      t2=p->w/t2;
    }
  p->ap=p->w/t*p->w/t/p->al/p->tfin*t2;
  p->bp=(p->w-p->dp)*(p->w-p->dp)/p->al/p->tfin;

  /* Generation time function */
  j=inputint(qstring,"pattype=", 0);
  p->gtime=(j==1)?gtliss:((j==2)?gtdisk:gtnorm);
  j=inputint(qstring,"xwave=", 1);
  p->xpat=(j==0)?sino:((j==1)?triang:sawt);
  j=inputint(qstring,"ywave=", 1);
  p->ypat=(j==0)?sino:((j==1)?triang:sawt);

  /* Set temporary file names */
  sprintf(gfile,"gen%d.ps",getpid()); p->gout=gfile;
  sprintf(hfile,"his%d.ps",getpid()); p->gtou=hfile;
  sprintf(sfile,"surf%d.ps",getpid()); p->sout=sfile;
  sprintf(lfile,"last%d.ps",getpid()); p->lout=lfile;
  sprintf(ffile,"flux%d.ps",getpid()); p->fout=ffile;
  if (p->gtime!=gtnorm)
    {
      sprintf(pfile,"patt%d.ps",getpid());
      p->pout=pfile;
    }
  else
    p->pout=NULL;
  if (p->cyc)
    {
      sprintf(tfile,"temp%d.ps",getpid()); p->tout=tfile;
    }
  else
    p->tout=NULL;

  return NULL;
}


/*++++++++++++++++++++++++++++++++++++++
  Search through template file for strang, printing as we search.

  FILE *infile The template file.

  char *strang String to search for.
  ++++++++++++++++++++++++++++++++++++++*/

void printto (FILE *infile, char *strang)
{
  char *a, temp[150];
  do
    {
      a = fgets (temp, 150, infile);
      if (a && mystrcmp (temp, strang))
        printf ("%s", temp);
    }
  while (a && mystrcmp (temp, strang));
}


/*++++++++++++++++++++++++++++++++++++++
  Search through file for strang.

  FILE *infile File to search through.

  char *strang String to search for.
  ++++++++++++++++++++++++++++++++++++++*/

void skipto (FILE *infile, char *strang)
{
  char *a, temp[150];
  do
    {
      a = fgets (temp, 150, infile);
    }
  while (a && mystrcmp (temp, strang));
}


/*++++++++++++++++++++++++++++++++++++++
  Generate a single postscript graph of various types.

  struct param *p EBsurf parameter structure.

  double *X The array of z-values.

  double *t The array of time step times.

  double *T The big array of temperatures.

  int msg Which graph to generate.
  ++++++++++++++++++++++++++++++++++++++*/

void output (struct param *p, double *X, double *t, double *T, int msg)
{
  int i,j,cyc;
  double tmp,maxi(),mini(),bot,top,time,scale;
  FILE *outs,*inbo;
  int steps,plots,flags;
  char bigs[81],inf[80],fname[80]=OUTPUTDIR;

  flags=((msg==TMPS||msg==ONEC)?1:0)|        /* 1: 3D */
    ((msg==SURF||msg==TMPS||msg==ONEC)?2:0)| /* 2: y or z axis "temperature" */
    ((msg==SURF||msg==FLUX||msg==GENT)?4:0); /* 4: first axis is time */
  strcpy(inf,PSDIR);
  if(flags&1)
    strcat (inf,"onedim3.ps");
  else
    strcat (inf,"onedim.ps");

  switch(msg)
    {
    case GENZ:
      strcpy(bigs,"Beam heat generation <I>G(z)</I>");
      strcat(fname,p->gout);
      break;
    case GENT: 
      strcpy(bigs,"Normalized beam history <I>G(t)</I>");
      strcat(fname,p->gtou);
      break;
    case SURF:
      strcpy(bigs,"Surface temperatures");
      strcat(fname,p->sout);
      break;
    case TMPS:
      strcpy(bigs,"All temperatures");
      strcat(fname,p->tout);
      break;
    case ONEC:
      strcpy(bigs,"Last cycle temperatures");
      strcat(fname,p->lout);
      break;
    case FLUX:
      strcpy(bigs,"Bottom heat flux");
      strcat(fname,p->fout);
      break;
    case PATT:
      strcpy(bigs,"Beam pattern");
      strcat(fname,p->pout);
      break;
    }

  printf("<TD><A HREF=\"%s%s\">%s</A></TD>\n", OUTPUTURLDIR,
         fname+strlen(OUTPUTDIR), bigs);
  if (!(outs=fopen(fname,"w")))
    exit(printf("No output file.\n"));
  if (!(inbo=fopen(inf,"r")))
    exit(printf("No input file.\n"));

  for (i=0;i<19;i++)
    {
      if (!(fgets(bigs,80,inbo)))
        exit(printf("Input error\n"));
      fputs(bigs,outs);
    }
  cyc = (p->cyc)?p->cyc:1;
  fflush (outs);

  steps = (flags&4) ? p->ts : p->zs; /* steps in each plot, plots=# of plots */
  plots = ((flags&4)&&msg!=GENT) ? cyc :
    1 + p->ts*((msg==TMPS) ? cyc : ((msg==GENZ||msg==GENT) ? 0 : 1));
  scale = (msg==FLUX) ? p->k/(p->w-X[p->zs-1])/1.e7 :
    ((msg==GENZ) ? p->k/p->al/1.e10 : 1);

  top = scale * ((msg==FLUX) ? (maxi(T+p->zs-1,p->ts*cyc+1,p->zs+1)-p->Tin) :
                 ((msg==GENT) ? 1 : maxi(T,(steps+1)*plots,1)));
  bot = (flags&2) ? mini(T,top,(steps+1)*plots,1) : 0;/* Minimum, maximum of */
                                                      /* abscissa */
  time = p->tfin * ((msg==GENT||msg==ONEC) ? 1 : cyc);
  if (flags&2)
    {
      top-=273.15; bot-=273.15;
    }
  if (msg==PATT)
    {
      time=p->plen; bot=-.5*p->pwid; top=.5*p->pwid;
    }

  if (flags&1)
    {
      fprintf (outs,"/lims { 0 %lf 0 %lf %lf %lf } def\n", time, p->w, bot,
               top-bot);
      fprintf (outs,"/labels { (t, sec) (z, cm) (Temperature, C) } def\n");
    }
  else
    {
      fprintf (outs,"/lims { %lf %lf %le %le } def\n",
              (msg==PATT)?-.5*time:0.,(msg==GENZ)?p->w:time,bot,top-bot);
      fprintf (outs,"/labels { (%s) (%s) } def\n",
               (msg==GENZ)?"depth, cm":"time, s",
               (msg==GENZ)?"Heat generation, kW/cm^2":
               ((msg==GENT)?"Normalized heat flux":
                ((msg==FLUX)?"Heat flux, W/cm^2":"Temperature, C")));
    }
  fflush (outs);

  for (i=0; i < ((flags&1) ? 136 : 115); i++)
    {
      if (!(fgets(bigs,80,inbo)))
        exit (printf ("Input read error...\n"));
      fputs(bigs,outs);
    }

  if (msg==PATT)
    { /* Create pattern graph */
      if(p->gtime==gtdisk)
        { /* Disk pattern */
          fprintf(outs,"%lf %lf gom\n",.5*p->plen,0);
          for(i=1;i<=p->ts;i++)
            fprintf(outs,"%lf %lf gol\n",
                    .25*p->plen*(cos(2.*M_PI*t[i]*p->xfrq) +
                                 cos(2.*M_PI*t[i]*p->yfrq)),