Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Posting-Version: notesfiles
Path: utzoo!watmath!clyde!bonnie!akgua!sdcsvax!sdcrdcf!hplabs!hp-pcd!hpcnoe!barrett
From: barrett@hpcnoe.UUCP (barrett)
Newsgroups: net.aviation
Subject: Re: Orphaned Response
Message-ID: <1200005@hpcnoe.UUCP>
Date: Tue, 16-Apr-85 19:53:00 EST
Article-I.D.: hpcnoe.1200005
Posted: Tue Apr 16 19:53:00 1985
Date-Received: Sun, 7-Apr-85 05:07:52 EST
References: <-26300@mss.UUCP>
Organization: Hewlett-Packard - Fort Collins, CO
Lines: 105
Nf-ID: #R:mss:-26300:hpcnoe:1200005:37777777600:3506
Nf-From: hpcnoe!barrett    Apr  2 16:53:00 1985

/*
 * Here is a program to play with density altitude:
 *
 * I have been working on this from time to time and just recently coded
 * this in C from my HP 15-C calculator program.  I havn't looked at it
 * for awhile, so I don't remember if it works or not.  The formula used
 * are derived from the universal gas law pv=nRT, and change in pressure
 * with altitude.  I have the derivations on paper, and translation to 
 * ascii is painful enough that I have never done so.  From what I can 
 * find out, the results should be within 1% accuracy (in an ideal universe)
 * up to about 38000 ft where the standard lapse rate no longer applies.
 *
 * Dave Barett
 * hplabs!hp-dcd!barrett		April 1, 1985
 */

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

#define AbsZero      (-273.16)
#define cent(Kelvin) ()
#define kelvin(Cent) ((Cent) - AbsZero)		      /* Centigrade to kelvin */
#define meters(feet) ((feet)*12*2.54/100) 	 /* meters to feet conversion */
#define FtoC(F)      (5/9*((F) - 32))
#define CtoF(C)      ((C) * 9/5 + 32)

#define LapseRate (-0.0065)  /* C/m */        /* Standard free air lapse rate */
#define T0        288.16     /* K */     /* Standard Temperature at sea level */
#define g         (-9.80655) /* m / s^2 */     /* Earth's gravitational const */
#define R         8.31432    /* kg m^2 / s^2*mol*K */  /* Universal Gas Const */
#define M         0.028966   /* kg / mol */            /* Molar weight of air */
#define P0        1013.25    /* mb= kg / 100*m*s^2 *//* pressure at sea level */
#define y0        0          /* m */         	/* true altitude at sea lavel */

#define StdTemp(y) (T0 + LapseRate*(y))  /* K */      /* Standard temperature */
#define StdAlt(T)  (((T)-T0)/LapseRate)  /* m */          /* Inverse of above */

#define Q         (g*M/(LapseRate*R))			     /* dimensionless */
#define k         (P0/pow(StdTemp(y0), Q))          /* mb / K */ /* d P / dy */
#define mu	  (100*k*M/R)	    	  /* kg / m^3 */        /* d rho / dy */


main()
{
    double alt, temp, press;
    double DensityAltitude(), stddensity(), stdpressure(), PressureAltitude();

    /*
    printf("y0   = %.6f\n", y0);
    printf("b    = %.6f\n", LapseRate);
    printf("a    = %.6f\n", T0);
    printf("g    = %.6f\n", g);
    printf("M    = %.6f\n", M);
    printf("R    = %.6f\n", R);
    printf("P0   = %.6f\n", P0);
    printf("y0   = %.6f\n", y0);
    printf("Q    = %.6f\n", Q);
    printf("k    = %.6e\n", k);
    printf("mu   = %.6e\n", mu);
    */

    printf(" Input alt(ft), Temp(C)?");
    scanf("%lf %lf", &alt, &temp);
    alt = meters(alt);
    printf(" alt          = %.6f m\n", alt);
    printf(" temp         = %.6f C\n", temp);
    printf(" Std Temp     = %.2f C\n", cent(StdTemp(alt)));
    printf(" Density      = %.4f kg / m^3\n", stddensity(alt));
    printf(" Pressure     = %.2f mb\n", stdpressure(alt, temp));
    printf(" Density alt  = %.2f ft\n", feet(DensityAltitude(alt, temp)));
}

double stddensity(altitude) 
   double altitude;
{
   return(mu*pow(StdTemp(altitude), Q-1));
}
   
double stddensityalt(density)
   double density;
{
   return(StdAlt(pow(density/mu, 1/(Q-1))));
}

double stdpressure(altitude) 
   double altitude;
{
   return(k*pow(StdTemp(altitude), Q));
}
   
double stdpressalt(pressure)
   double pressure;
{
   return(StdAlt(pow(pressure/k, 1/Q)));
}

double DensityAltitude(altitude, temperature)
   double altitude, temperature;
{
   return(
      stddensityalt(
	 M/R * 100 * stdpressure(altitude) / kelvin(temperature)
      )
   );
}