standardatmosphere.m

matlab程序

                 "Altitude, Km"},
    PlotRange->All, 
    Axes->False]

(* values from the 'U.S. Standard Atmosphere, 1976' model *)
$Density=(* (Kilo Gram/Meter^3) *){1.9311, 1.8491, 1.7697, 1.6930, 1.6189, 
    1.5473, 1.4782, 1.4114, 1.3470, 1.2849, 1.2250, 1.1673, 1.1117, 1.0581, 
    1.0066, 9.5695 10^-1, 9.0925 10^-1, 8.6340 10^-1, 8.1935 10^-1, 
    7.7704 10^-1, 7.3643 10^-1, 6.9747 10^-1, 6.6011 10^-1, 6.2431 10^-1, 
    5.9002 10^-1, 5.5719 10^-1, 5.2579 10^-1, 4.9576 10^-1, 4.6706 10^-1, 
    4.3966 10^-1, 4.1351 10^-1, 3.8857 10^-1, 3.6480 10^-1, 3.3743 10^-1, 
    3.1194 10^-1, 2.8838 10^-1, 2.6660 10^-1, 2.4646 10^-1, 2.2786 10^-1, 
    2.1066 10^-1, 1.9476 10^-1, 1.6647 10^-1, 1.4230 10^-1, 1.2165 10^-1, 
    1.0400 10^-1, 8.8910 10^-2,7.5715 10^-2, 6.4510 10^-2, 5.5006 10^-2, 
    4.6938 10^-2, 4.0084 10^-2, 3.4257 10^-2, 2.9298 10^-2, 2.5076 10^-2, 
    2.1478 10^-2, 1.8410 10^-2, 1.5792 10^-2, 1.3555 10^-2, 1.1573 10^-2, 
    9.8874 10^-3, 8.4634   10^-3, 7.2579 10^-3, 5.3666 10^-3, 3.9957 10^-3, 
    2.9948 10^-3, 2.2589 10^-3, 1.7142 10^-3, 1.3167 10^-3, 1.0269 10^-3, 
    8.0562 10^-4, 6.3901 10^-4, 5.0445 10^-4, 3.9627 10^-4, 3.0968 10^-4, 
    1.6321 10^-4, 8.2829 10^-5, 3.9921 10^-5, 1.8458 10^-5, 8.2196 10^-6, 
    7.5641 10^-6,6.958 10^-6, 3.416 10^-6, 1.393 10^-6, 5.604 10^-7, 
    9.708 10^-8,   2.222 10^-8, 8.152 10^-9, 3.831 10^-9, 2.076 10^-9, 
    1.233 10^-9, 7.815 10^-10, 5.194 10^-10, 3.581 10^-10, 2.541 10^-10, 
    1.846 10^-10, 1.367 10^-10, 7.858 10^-11, 4.742 10^-11, 2.971 10^-11, 
    1.916 10^-11, 1.264 10^-11, 8.503 10^-12, 5.805 10^-12, 4.013 10^-12, 
    2.803 10^-12, 1.184 10^-12, 5.215 10^-13, 2.384 10^-13, 1.137 10^-13, 
    5.712 10^-14, 3.070 10^-14, 1.788 10^-14, 1.136 10^-14, 7.824 10^-15, 
    5.759 10^-15, 4.453 10^-15, 3.561 10^-15 };

AtmosphericPlot[MeanDensity, opts___?OptionQ]:=
  ListPlot[
    Transpose[{Log[10, $Density], $Alts/1000}],
    opts,
    PlotJoined->True,
    Frame->True,
    FrameLabel->{"Log[10, MeanDensity]",
                 "Altitude, Km"},
    RotateLabel->True, 
    PlotRange->All,
    Axes->False]

$SoundSpeed=(*(Meter/Second) *) {358.99, 357.16, 355.32, 353.48, 351.63, 
    349.76, 347.89, 346.00, 344.11, 342.21, 340.29, 338.37, 336.43, 334.49, 
    332.53, 330.56, 328.58, 326.59, 324.59, 322.57, 320.55, 318.50, 316.45, 
    314.39, 312.31, 310.21, 308.11, 305.98, 303.85, 301.7, 299.53, 297.35, 
    295.15, 295.07, 295.07, 295.07, 295.07, 295.07, 295.07, 295.07, 295.07, 
    295.07, 295.07, 295.07, 295.07, 295.07, 295.07, 296.38, 297.05, 297.72, 
    298.39, 299.06, 299.72, 300.39, 301.05, 301.71, 302.37, 303.02, 304.67, 
    306.49, 308.30, 310.10,313.67, 317.19, 320.67, 324.12, 327.52, 329.80, 
    329.80, 328.81, 325.43, 322.01, 318.56, 315.07, 306.19, 297.06, 289.40, 
    282.54, 275.52, 274.81} ;

AtmosphericPlot[SoundSpeed, opts___?OptionQ]:=
  ListPlot[
    Transpose[{$SoundSpeed, Take[$Alts/1000, Length[$SoundSpeed]]}],
    opts,
    PlotJoined->True,
    Frame->True,
    FrameLabel->{"SoundSpeed, m/s",
                 "Altitude, Km"},
    PlotRange->All,
    Axes->False];

$DynamicViscosity=(* (Newton Second/Meter^2) *) {19.422, 19.273, 19.123, 
      18.972, 18.820, 18.668, 18.515, 18.361, 18.206, 18.050, 17.894, 17.737, 
      17.579, 17.420, 17.260, 17.099, 16.938, 16.775, 16.612, 16.448, 16.282, 
      16.116, 15.949, 15.781, 15.612, 15.442, 15.271, 15.099, 14.926, 14.752, 
      14.577, 14.400, 14.223, 14.216, 14.216, 14.216, 14.216, 14.216, 14.216,
      14.216, 14.216, 14.216, 14.216, 14.216, 14.216, 14.216, 14.267, 
      14.322, 14.376, 14.430, 14.484, 14.538, 14.592, 14.646, 14.699, 14.753, 
      14.806, 14.859, 14.992, 15.140, 15.287, 15.433, 15.723, 16.009, 16.293, 
      16.573, 16.851, 17.037, 17.037, 16.956, 16.600, 16.402, 16.121, 15.837, 
      15.116, 14.377, 13.759, 13.208, 12.647, 12.590} 10^-6;

AtmosphericPlot[DynamicViscosity, opts___?OptionQ]:=
  ListPlot[
    Transpose[{10^6 $DynamicViscosity,
	 Take[$Alts/1000, Length[$DynamicViscosity]]}],
    opts,
    PlotJoined->True,
    Frame->True,
    FrameLabel->{"DynamicViscosity, 10^-6 N s/m^2",
                 "Altitude, Km"},
    PlotRange->All,
    Axes->False];

$KinematicViscosity=(* ((Meter^2)/Second) *){1.0058 10^-5, 1.0423 10^-5, 
    1.0806 10^-5, 1.1206 10^-5, 1.1225 10^-5, 1.2065 10^-5, 1.2525 10^-5, 
    1.3009 10^-5, 1.3516 10^-5, 1.4048 10^-5, 1.4607 10^-5, 1.5195 10^-5, 
    1.5813 10^-5, 1.6463 10^-5, 1.7147 10^-5, 1.7868 10^-5, 1.8628 10^-5, 
    1.9429 10^-5, 2.0275 10^-5, 2.1167 10^-5, 2.2110 10^-5, 2.3107 10^-5, 
    2.4161 10^-5, 2.5278 10^-5, 2.6461 10^-5, 2.7714 10^-5, 2.9044 10^-5, 
    3.0457 10^-5, 3.1957 10^-5, 3.3553 10^-5, 3.5251 10^-5, 3.7060 10^-5, 
    3.8988 10^-5, 4.2131 10^-5, 4.5574 10^-5, 4.9297 10^-5, 5.3325 10^-5, 
    5.7680 10^-5, 6.2391 10^-5, 6.7485 10^-5, 7.2995 10^-5, 8.5397 10^-5, 
    9.9901 10^-5, 1.1686 10^-4,1.3670 10^-4, 1.5989 10^-4, 1.8843 10^-4, 
    2.2201 10^-4, 2.6135 10^-4, 3.0743 10^-4, 3.6135 10^-4, 4.2439 10^-4, 
    4.9805 10^-4, 5.8405 10^-4, 6.8437 10^-4, 8.0134 10^-4, 9.3759 10^-4, 
    1.0962 10^-3, 1.2955 10^-3, 1.5312 10^-3, 1.8062 10^-3, 2.1264 10^-3, 
    2.9297 10^-3, 4.0066 10^-3, 5.4404 10^-3, 7.3371 10^-3, 9.8305 10^-3, 
    1.2939 10^-2, 1.6591 10^-2, 2.1047  10^-2, 2.6104 10^-2, 3.2514 10^-2, 
    4.0682 10^-2, 5.1141 10^-2, 9.2617 10^-2, 1.7357 10^-1, 3.4465 10^-1, 
    7.1557 10^-1, 1.5386, 1.6645}  ;

AtmosphericPlot[KinematicViscosity, opts___?OptionQ]:=
  ListPlot[
    Transpose[{Log[10, $KinematicViscosity],
                   Take[$Alts/1000, Length[$KinematicViscosity]]}],
    opts,
    PlotJoined->True,
    Frame->True,
    FrameLabel->{"Log[10, KinematicViscosity]",
                 "Altitude, Km"},
   PlotRange->All,
   Axes->False]

$ThermalConductivity=(* (Joule/(Meter Second Kelvin)) *){2.7882, 2.7634, 
    2.7384, 2.7134, 2.6884, 2.6632, 2.6380, 2.6126, 2.5872, 2.5618, 2.5326, 
    2.5106, 2.4849, 2.4591, 2.4333, 2.4073, 2.3813, 2.3552, 2.3290, 2.3028, 
    2.2765, 2.2500, 2.2236, 2.1970, 2.1703, 2.1436, 2.1168, 2.0899, 2.0630, 
    2.0359, 2.0088, 1.9816, 1.9543, 1.9533, 1.9533, 1.9533, 1.9533,1.9533, 
    1.9533, 1.9533, 1.9533, 1.9533, 1.9533, 1.9533, 1.9533, 1.9533, 1.9611, 
    1.9695, 1.9778, 1.9862, 1.9945, 2.0029, 2.0112, 2.0195, 2.0278, 2.0361, 
    2.0443, 2.0526, 2.0733, 2.0963, 2.1193, 2.1422, 2.1878, 2.2331, 2.2781, 
    2.3229, 2.3764, 2.3973, 2.3973, 2.3843, 2.3400, 2.2955, 2.2508, 2.2058, 
    2.0926, 1.9780, 1.8834, 1.8001, 1.7126, 1.7078};

AtmosphericPlot[ThermalConductivityCoefficient, opts___?OptionQ]:=
 ListPlot[
    Transpose[{$ThermalConductivity,
	 Take[$Alts/1000, Length[$ThermalConductivity]]}],
    opts,
    PlotJoined->True,
    Frame->True,
    FrameLabel->{"ThermalConductivityCoefficient,J/(m s K)",
                 "Altitude, Km"},
   PlotRange->All,
   Axes->False]

Scan[
 (
  MessageName[#, "cnvrt"] =
	"Unable to convert altitude `1` to a numerical multiple of Meter.";
  MessageName[#, "rng"] =
	"The altitude `1` falls outside of the 
	supported range from `2` to `3` Meters."
 )&,
	{CollisionFrequency, DynamicViscosity, GravityAcceleration, 
   KinematicViscosity, KineticTemperature, MeanDensity, MeanFreePath, 
   MeanParticleSpeed, MeanMolecularWeight, NumberDensity, Pressure, 
   PressureScaleHeight, SoundSpeed, ThermalConductivityCoefficient}
]

atmosphericDataInterpolation[data_, altitude_, unit_, object_] :=
   Module[{n, alts, low, high, pos, a, p, a1, a2, d1, d2,
	 alt = Convert[altitude, Meter]/Meter},
     If[!NumberQ[N[alt]],
	 Message[MessageName[object, "cnvrt"], altitude];
	 Return[$Failed]];
     n = Length[data];
     alts = Take[$Alts, n];
     {low, high} = {First[alts], Last[alts]};
     If[alt < low || alt > high,
	Message[MessageName[object, "rng"], altitude, low, high];
	Return[$Failed]
     ];
     pos = Position[alts, a_ /; a>=alt, 1, 1];
     If[MatchQ[pos, {{1}}] && alt==low,
        Return[First[data]*unit]];
     p = First[First[pos]];
     {a1, a2} = alts[[{p-1, p}]];
     {d1, d2} = data[[{p-1, p}]];
     (d1 + (alt-a1)/(a2-a1) * (d2-d1))*unit
		]

End[];

EndPackage[];

(* :Examples:
	error examples:
		CollisionFrequency[1]
		MeanParticleSpeed[-6000 Meter]
		GravityAcceleration[1000001 Meter]		
		DynamicViscosity[85501 Meter]
*)

(* :Examples:
	AtmosphericPlot[CollisionFrequency]
	AtmosphericPlot[DynamicViscosity]
		(* DynamicViscosity does not decrease monotonically
			with increasing altitude. *)
	AtmosphericPlot[GravityAcceleration]
	AtmosphericPlot[KinematicViscosity]
	AtmosphericPlot[KineticTemperature]
		(* KineticTemperature does not increase monotonically
                        with increasing altitude. *)
	AtmosphericPlot[MeanDensity]
	AtmosphericPlot[MeanFreePath]
	AtmosphericPlot[MeanMolecularWeight]
	AtmosphericPlot[MeanParticleSpeed]
		(* MeanParticleSpeed does not increase monotonically
                        with increasing altitude. *)
	AtmosphericPlot[NumberDensity]
	AtmosphericPlot[Pressure]
	AtmosphericPlot[PressureScaleHeight]
		(* PressureScaleHeight does not increase monotonically
                        with increasing altitude. *)
	AtmosphericPlot[SoundSpeed]
		(* SoundSpeed does not decrease monotonically 
			with increasing altitude. *)
	AtmosphericPlot[ThermalConductivityCoefficient]
		(* ThermalConductivityCoefficient does not decrease
                        monotonically with increasing altitude. *)
*)

(* :Examples:
	Plot[ Evaluate[CollisionFrequency[x Meter] Second],
			{x, -5000, 1000000}]
	Plot[ Evaluate[DynamicViscosity[x Meter] Meter^2/(Newton Second)],
			{x, -5000, 85500}]
	Plot[ Evaluate[GravityAcceleration[x Meter] Second^2/Meter],
			{x, -5000, 1000000}]
	Plot[ Evaluate[KinematicViscosity[x Meter] Second/Meter^2],
			{x, -5000, 85500}]
	Plot[ Evaluate[KineticTemperature[x Meter] 1/Kelvin],
			{x, -5000, 1000000}]
	Plot[ Evaluate[MeanDensity[x Meter] Meter^3/Kilogram],
			{x, -5000, 1000000}]
	Plot[ Evaluate[MeanFreePath[x Meter] 1/Meter],
			{x, -5000, 1000000}]
	Plot[ Evaluate[MeanMolecularWeight[x Meter] Kilo Mole/Kilogram],
			{x, -5000, 1000000}]
	Plot[ Evaluate[MeanParticleSpeed[x Meter] Second/Meter], 
        		{x, -5000, 1000000},
	   AxesLabel -> {"m", "m/sec"}]
	Plot[ Evaluate[NumberDensity[x Meter] Meter^3],
			{x, -5000, 1000000}]
	Plot[ Evaluate[Pressure[x Meter] 1/(Milli Bar)],
			{x, -5000, 1000000}]
	Plot[ Evaluate[PressureScaleHeight[x Meter] 1/Meter],
			{x, -5000, 1000000}]
	Plot[ Evaluate[SoundSpeed[x Meter] Second/Meter],
			{x, -5000, 85500}]
	Plot[ Evaluate[ThermalConductivityCoefficient[x Meter]
			 (Meter Second Kelvin)/Joule],
			{x, -5000, 85500}]
*)