/*-------------------------------------------------------------------\
|  Thermal Index Calculator (ti)                       Version 1.60  |
|  written by Kevin Ford                               1-28-94       |
|  version 1 written Feb. 1992                                       |
|								     |
|  VAX/VMS Compatibility added 06/30/93 by Kim Reniska               |
|                                                                    |
|  Six modes of operation:				             |
|  Interactive:   						     |
|     ti                read from .uac file (text)                   |
|     ti -1		read from .uac file			     |
|     ti -3	        read from archived file			     |
|  Batch							     |
|     ti    filename station hightemp maxz  ; read from .uac (text)  |
|     ti -1 filename station hightemp maxz  ;  read from .uac file   |
|     ti -3 filename date hightemp maxz  ;  read from archive file   |
\-------------------------------------------------------------------*/

/* must use cc ti.c -lm to compile */

#include <stdio.h>
#include <math.h>
#include <ctype.h>
#include <string.h>
#ifdef VMS
#define O_RDONLY 0
#include <unixio.h>
#else
#include <fcntl.h>
#endif

#define t0  (273.15)      /* convert from Kelvin to Celcius */
#define RG (29.286)      
#define GD (-0.00974)      /* dry adiabatic lapse rate, deg. C per meter */
#define G0 (-0.0065)       /* standard lapse rate */
#define CFM  (0.3048)        /* convert from feet to meters */
#define spdcvt (1.9458)   /* convert from m/s to knots */

/* stored data structure for .uac files */

typedef struct { 
  char ident[12], loc[12], loc2[12];
  long idcode;
  float lat,lon;
  long ml,sl,wl;   /* num. of mand. levels, sig. levels, wind levels */
  float mp[20], mh[20], mt[20], mdp[20], mdir[20], mspd[20];
  float sp[50], st[50], sdp[50];
  float wp[50], wdir[50], wspd[50];
} filerec;
filerec rec;

/* raw data from the file, for each level */

float p[100],           /* pressure in millibars */
      t[100],           /* temperature in degrees Kelvin */
      dp[100];          /* dewpoint in degrees Kelvin    */
int   ddir[100],        /* wind direction, degrees */
      dspd[100];        /* wind speed in knots */

/* auxilliary 1000 foot interval wind data from the file */

int winddir[160],       /* wind direction, degrees */
    windspd[160];       /* wind speed, knots       */

/* calculated data for each pressure level */

float h[100],           /* height of level in meters MSL  */
      tv[100];          /* virtual temperature in degrees Celcius  */

/* These arrays contain calculated data at 500 foot intervals */

float ti[300],          /* thermal index     */
      alt[300];         /* altimeter reading */
int   dir[300],         /* wind direction, degrees */
      spd[300];         /* wind speed in knots */

int levels;             /* total number of pressure levels */

/* calculate the virtual temp given the temp, dewpoint and presure */

float vtemp(t,dp,p)
  float t,dp,p;
{ float e;
  e=6.11 * exp(5423.0*(1.0/t0 - 1.0/dp));
  return(t/(1.0-.378*e/p));
}

/* altitude for a given pressure in the standard atmosphere */

float stdz(p) 
  float p;
{
  return((-288.0/G0) * (1.0- pow(p/1013.25,-RG*G0)));
}

/* convert a string to an integer */

int strtoint(buf)
  char *buf;
{
  int d; char w;
  d=0;
  while ((w= *buf++)>0) d=10*d+w-'0';
  return(d);
}

/* decode a coded temperature/dewpoint from an text data file */

tempcvt(n,t,dp)
  int n;         /* coded form */
  float *t, *dp; /* decoded temp & dewpoint */
{
  int q; float h;
  q=n/100;
  h=((float) q)/10.0;
  if (q%2==0) *t=t0+h; else *t=t0-h;
  q=n%100;
  h=(float) q;
  if (q<56) *dp = *t - h/10.0;
    else if (q<99) *dp = *t - h+50.0;
      else *dp = t0 - 99.9;
}

/* virtual temperature at altitude a */

float tempalt(a)
  float a;
{
  int i;
  if (a<=h[0]) return(tv[0]);
  if (a >= h[levels-1]) return(tv[levels-1]);
  for (i=1; h[i]<a; i++);
  return(tv[i-1]+(tv[i]-tv[i-1])*(a-h[i-1])/(h[i]-h[i-1]));
}

main(argc, argv)
  int argc;
  char *argv[];
 {
  int index,itemp,fscanf(),found,pr,p0,n,m,WND,hi,lo,ht,hv,i,j,
      fd, nid, stat, ti1, ti3, batch, trig[200], cb;
/*  long lseek(), read(), open();  FILE *fopen(),*fp; */
  FILE *fp;
  char buf[200],fname[100], id[10], date[20], date2[20],
     block[300][10], station[20];
  float salt,statalt,ftemp,htt,hd,k,ta,z,a,maxz;

/*--------------------- Process command line arguments ----------------*/

  ti1 = 0; ti3=0; batch=0; i=1;
  if (argc>1) {
     if (strcmp(argv[1],"-3")==0) { ti3=1; i++; }
     if (strcmp(argv[1],"-1")==0) { ti1=1; i++; }
  }
  if (argc>i) {
    if (argc<i+4) { printf("ti error: Too few arguments\n"); exit(-1); }
    strncpy(fname, argv[i], 99); fname[99]=0; i++;
    if (ti3) { strncpy(date,argv[i],10); date[10]=0; }
    else { strncpy(station,argv[i],6); station[6]=0; }
    i++;
    sscanf(argv[i],"%d",&ht); i++;
    sscanf(argv[i],"%f",&maxz);
    batch=1;
  }

/*------------ Print headers and open file ------------- */
  
  if (!batch) {
    printf("Thermal Index Calculator\n");
    printf("     by Kevin Ford\n\n");
    printf("Reading files in the wxp ");
    if (ti3) printf("text archive format (one station per file)\n");
    if (ti1+ti3==0) printf("text .uac format (one time per file)\n");
    if (ti1) printf(".uac format (non-text)\n");
      printf("\nEnter filename : ");
      scanf("%s",fname);
  }
  if (!ti1) {
    if ((fp=fopen(fname,"r")) == NULL) {
        printf("ti error: cannot open file %s.\n",fname); exit(-1);
    }
  }
  else if ((fd = open(fname, O_RDONLY,0)) == -1) {
     printf("ti error: cannot open file %s.\n",fname); exit(-1);
  }

/*--------- Get date/station and read data for text data file --------------*/

  if (!ti1) {
    if (!batch) {
       if (ti3) { printf("Enter date (YYMMDDHH) : "); scanf("%s",date); }
       else { printf("Enter station : "); scanf("%s",station); }
    }
    found=0;
    while (!found) {
      if (fscanf(fp,"%s",buf)==EOF) found=2;
      if ((strncmp(buf,"DATE",4)==0) && (ti3)) {
        fscanf(fp,"%s",date2);
        if (strncmp(date2,date,8)==0) found=1;
      }
      if ((!ti3) && (strncmp(buf,station,3) == 0)) found=1;
    }
    if (found!=1) {
      if (ti3) printf("ti error: data for %s not found.\n",date); 
      else printf("ti error: data for %s not found.\n",station);
      exit(-1);
    }
    if (ti3) fscanf(fp,"%s",station);
    for (i=0; block[i-1][0]!='$'; i++) fscanf(fp,"%s",block[i]);
    fclose(fp);
  }   /* if !ti1 */

/*--------- Get station id and read upper air data from .uac file -------*/

  else {
  lseek(fd,4284L,0); nid=0;
  while ((read(fd,&rec,sizeof(filerec))) == sizeof(filerec))
    strcpy(block[nid++],rec.ident);
  if (nid==0) {
     printf("ti error: file %s contains no sounding data.\n",fname);
     exit(-1);
  }

  if (!batch) {
    printf("Stations:\n\n");
    for (i=0; i<nid; i++) {
       block[i][4]=0;
       printf("%5s",block[i]); if (i%15==14) printf("\n");
    }
    for (stat= -1; stat<0;) {
      printf("\n\nEnter station id: "); scanf("%s",station);
      for (i=0; station[i]!=0; i++) 
        if ((station[i]>='a')&&(station[i]<='z')) station[i] -= ('a'-'A');
      for (j=0; j<nid; j++) if (strcmp(station,block[j])==0) { stat=j;break; }
    }
  }
  else {
    stat= -1;
    for (i=0; station[i]!=0; i++) 
      if ((station[i]>='a')&&(station[i]<='z')) station[i] -= ('a'-'A');
    for (j=0; j<nid; j++) if (strcmp(station,block[j])==0) { stat=j;break; }
    if (stat<0) {
      printf("ti error: station %s data not found.\n",station);
      exit(-1);
    } 
  }

  lseek(fd,4284L+stat*sizeof(filerec),0);
  read(fd,&rec,sizeof(filerec));
  close(fd);
  }  /* if ti1 */

/*------------ Read in station altitude from data file ------------------*/

  statalt = 0.0;
  if ((fp = fopen("stations","r")) != NULL) {
    while (fscanf(fp,"%s %f\n",id,&salt) != EOF) 
      if (strcmp(id,station)==0) { statalt = salt; break; }
    fclose(fp);
  }

/*--------------- Parse data from text data file ---------------------*/

  if (!ti1) {

              /* read in significant levels */

  p0= -999;
  levels=0; j=53;
  while ((block[j][0]!='X')&&(block[j][0]!='-')) {
    pr = strtoint(block[j]);  if (block[j][0]=='0') pr+=1000;
    if (j==53) p0=pr;  /* sfc pressure */
    p[levels] = (float) pr;
    dspd[levels] = -999; ddir[levels] = -999;
    j++;
    if (block[j][0]!='X') {
       tempcvt(strtoint(block[j]),&t[levels],&dp[levels]); j++;
    }
    else { t[levels]= -999.0; dp[levels] = -999.0; }
    levels++;
  }
  j++; while (block[j][0]=='X') j++;
  index=j;
  if (levels<1) {printf("ti error: sfc pressure missing.\n");exit(-1); }

            /* read in mandatory levels */
 
  WND=0;       
  for (j=0; j<53; j+=3) {
    if (block[j][0]!='X') {
      n = strtoint(block[j]);
      pr = n/1000; pr*=10; if (pr==0) pr=1000;
      if (pr==880) pr=n%1000;                   /* level of tropopause */
      else if (pr==770) { pr=n%1000; WND=1; }   /* level of max. wind  */
      if (pr<p0) {
        p[levels] = (float) pr;
        if (!WND) {
          if (block[j+1][0]!='X') 
	    tempcvt(strtoint(block[j+1]),&t[levels],&dp[levels]);
          else { t[levels]= -999.0; dp[levels] = -999.0; }
        }
	else {t[levels]= -999.0; dp[levels] = -999.0; }
        if (block[j+2-WND][0]!='X') { 
          n=strtoint(block[j+2-WND]); 
          ddir[levels] = (n/500)*5;
          dspd[levels] = n%500;
        }
        else {
          ddir[levels] = -999; dspd[levels] = -999;
        }
        levels++;
      }
    }
  }        

               /* read in winds aloft data */

  j=index;
  for (i=0; i<=150; i++) { windspd[i] = -999; winddir[i] = -999; }
  while (block[j][0]!='$') {
    n=strtoint(block[j]);
    if ((block[j+1][0]!='X') && (block[j+1][0]!='$')) {
      m=strtoint(block[j+1]);
      windspd[n] = m%500;
      winddir[n] = (m/500)*5;
    }
    j += 2;
    if (block[j-1][0]=='$') j--;
  }

  } /* if !ti1 */

/*------------------ Organize data from .uac file --------------------*/

  else {

  for (i=0;i<rec.sl;i++) {
    p[i]=rec.sp[i];
    t[i]=rec.st[i] + t0;
    dp[i]=rec.sdp[i] + t0;
    if (dp[i] < 100.0) dp[i] = t0 - 99.9;
    dspd[i]= -999;
    ddir[i]= -999;
  }
  for (i=0,j=rec.sl;i<rec.ml;i++,j++) {
    p[j]=rec.mp[i];
    t[j]=rec.mt[i] + t0;
    dp[j]=rec.mdp[i] + t0;
    if (dp[j] < 100.0) dp[j] = t0 - 99.9;
    dspd[j]= (int) (rec.mspd[i] * spdcvt + 0.5);
    ddir[j]= ((int) ((rec.mdir[i] + 2.5)/5.0)) * 5;
  }
  for (i=0;i<100;i++) { windspd[i]= -999.0; winddir[i]= -999.0; }
  for (i=0; i<rec.wl; i++) {
    j = (int) ((rec.wp[i]/CFM + 500.0)/1000.0);
    windspd[j] = (int) (rec.wspd[i] * spdcvt + 0.5);
    winddir[j] = ((int) ((rec.wdir[i] + 2.5)/5.0)) * 5;
  }
  levels=rec.ml+rec.sl;
  }  /* if ti1 */

/*-------------------- Sort data by pressure ------------------------*/

/* printf("raw data:\n");
   for (i=0; i<levels;i++) printf("%8.1f%8.1f%8.1f\n",p[i],t[i],dp[i]); 
*/

  for (i=0;i<levels-1;i++) for (j=i+1;j<levels;j++)
    if ((p[i]<p[j])||((p[i]==p[j])&&(dspd[i]>dspd[j]))) {
        ftemp=p[i];p[i]=p[j];p[j]=ftemp;
        ftemp=t[i];t[i]=t[j];t[j]=ftemp;
        ftemp=dp[i];dp[i]=dp[j];dp[j]=ftemp;
        itemp=dspd[i];dspd[i]=dspd[j];dspd[j]=itemp;
        itemp=ddir[i];ddir[i]=ddir[j];ddir[j]=itemp;
    }

  p[levels]=9999.0;
  for (i=0; i<levels; i++) 
    if ((t[i]<0.0)||(p[i]<9.9)||(p[i]-p[i+1]<.1)) {
      levels--;
      for (j=i; j<=levels; j++) {
        p[j]=p[j+1]; t[j]=t[j+1]; dp[j]=dp[j+1]; dspd[j]=dspd[j+1];
        ddir[j]=ddir[j+1];
      }
      i--;
    }

/* printf("\sorted data:\n");
   for (i=0; i<levels;i++) printf("%8.1f%8.1f%8.1f%8d%8d\n",
   p[i],t[i],dp[i], dspd[i],ddir[i]);  */

  if ((levels<2) || ((ti1) && (rec.mh[0] < -100.0))) { 
    printf("ti error: insufficient upper air data for ");
    if (ti3) printf("%s\n",date);
    else printf("%s\n",station);
    exit(-1);
  }

/*------------------ calculate heights of levels ------------------------*/

   /* get surface elevation */

  h[0]=statalt;

/*  printf("station elevation: %f\n",statalt);
  if (ti1) printf("lowest ml altitude: %f\n",rec.mh[0]);  */

  if (ti1) if ((h[0]-rec.mh[0])*(h[0]-rec.mh[0]) > 10.0) {
     printf("ti error: surface data missing\n");
     exit(-1);
  }

  for (i=0; i<levels;i++) tv[i] = vtemp(t[i],dp[i],p[i]);
  for (i=1; i<levels;i++) {
    ftemp = log(p[i-1]) - log(p[i]);
    if (tv[i]==tv[i-1])  h[i]=h[i-1]+RG*tv[i]*ftemp;
    else h[i]=h[i-1]+RG*(tv[i]-tv[i-1]) * ftemp/(log(tv[i])-log(tv[i-1]));
  }

/*--------------------- Get forecast hi, dp, max alt to calc TI ---------*/

    /* set dewpoint = average dewpoint in lowest 500 meters */

  z = 0.0; i=1; while (h[i] < h[0] + 500.0) {
    z += ((h[i]-h[i-1])*(dp[i]+dp[i-1])/2.0);
    i++;
  }
  if (i >= levels) hd=dp[0];
  else {
    ftemp = h[i] - h[i-1];
    a = dp[i-1]*(h[i]-(h[0]+500.0))/ftemp + dp[i]*(500.0+h[0]-h[i-1])/ftemp;
    z += (500.0+h[0]-h[i-1]) * (a + dp[i-1]) / 2.0;
    hd = z / 500.0;
  }

  if (!batch) {
    printf("\nEnter forecast high temp (F) : "); scanf("%d",&ht);
    printf("Enter maximum altitude to calculate TI : "); scanf("%f",&maxz);
  }
  ftemp = (float) (ht-32.0)/1.8+t0;
  if (ftemp<t[0]) { ftemp=t[0]; ht = (int) ((ftemp-t0)*1.8+32.5); }
  if (maxz*CFM<h[0]+305.0) maxz=(h[0]+305.0)/CFM;
  if (hd<200.0) hd = t[0]-5.0;
  if (hd>ftemp) hd=ftemp;
  cb = (int) ((ftemp - hd)*4.0 + .5) * 100;  /* est. cloud base */
  htt = vtemp(ftemp,hd,p[0]);

  printf("\nSoaring report from %s upper air data\n\n",station);
  printf("Forecast high : %3d F\n", ht);
  printf("Surface virtual temp = %4.1f C\n\n",htt-t0+.05);

  for (i=0; i<199; i++) {dir[i]= -999.0; spd[i]= -999.0; }

/*------------------ calculate thermal index ---------------------------*/

  lo=(int) (h[0]/CFM/500.0) + 1;
  hi=(int) ((maxz+250.0)/500.0);
  if (hi>120) hi=120;  /* don't go over 60,000 feet MSL */
  if (hi<lo+9) hi=lo+9;   /* do at least 5000 feet AGL */

  for (i=0; i<levels; i++)
    if ((ddir[i]>=0) && (dspd[i]>=0)) {
        j = (int) ((h[i]/CFM+250.0)/500.0);
        dir[j] = ddir[i]; spd[j] = dspd[i];
  }

  for (i=lo; i<=hi; i++) {
    k=500.0 * CFM * i;
    if ((i%2==0)&&(winddir[i/2]>=0)) {
      dir[i]=winddir[i/2]; spd[i]=windspd[i/2];
    }
    for (j=1; j<levels;j++) if ((h[j-1]<k)&&(k<=h[j])) {
        ta=tv[j-1] + (tv[j]-tv[j-1]) * (k-h[j-1])/ (h[j]-h[j-1]);
        ti[i] = ta - (htt + GD * (k-h[0]));
        trig[i] = ht + (int) ((ti[i]+3.0)*1.8+.5);
        break;
    }
  }

/*---------------------- print out tables ----------------------------*/

  printf("\n   P(mb)     H(ft)    Tv(C)    T(C)    DP(C)   ");
  printf("wind dir  wind spd\n\n");
  for (i=0; i<levels;i++) {
    printf(" %8.1f%9d%9.1f%9.1f",p[i],(int) (h[i]/CFM+.5), tv[i]-t0+.01,
       t[i]-t0+.01);
    if (dp[i]>180.0) printf("%9.1f",dp[i]-t0+.01);
        else printf("         ");
    if (dspd[i]>=0) printf("%9d%9d",ddir[i], dspd[i]);
    printf("\n");
  }
  if (p[0]<500.0) {
     printf("Insufficient upper air data for thermal index calculations.\n");
     exit(-1);
  }
  printf("\n\n\n  alt(ft MSL)     ti    wind dir   wind spd  trig temp (F)\n\n");
  for (i=lo; i<=hi; i++) {
    printf("%10d%10.1f%",i*500, ti[i]+.05);
    if ((spd[i]<0)||(dir[i]<0)) printf("                    ");
    else printf("%10d%10d",dir[i],spd[i]);
    printf(" %10d\n", trig[i]);
  }
  printf("\nEstimated cloud base: %d feet AGL\n",cb);
}