# Graphic manipulation - Psychrometric chart

I'm working on a psychrometric chart on Mathematica. The code is as follows:

(*psat is saturation pressure calculated by Hyland Wexler correlation*)
psat[T_?NumericQ] :=
Exp[(c8/(T + 273.15)) +
c9 + (c10*(T +
273.15)) + (c11*((T + 273.15)^2)) + (c12*((T +
273.15)^3)) + (c13*Log[T + 273.15])]
c8 = -5800.2206;
c9 = 1.3914993;
c10 = -0.04860239;
c11 = 0.41764768*10^-4;
c12 = -0.14452093*10^-7;
c13 = 6.5459673;
patm = 101325;(*Pa*)(*atm pressure*)
rn = 0.622;(*kg/kg*);(*molar ratio*)
cp = 1.004;(*KJ/kg K*);(*air specific heat*)
hvref = 2501;(*KJ/kg*);(*vapour enthalpy*)
cpda = 1.006;(*KJ/kg K*);(*specific heat dry air*)
cpwv = 1.84;(*KJ/kg K*);(*specific heat water vapour*)

(*Y is humidity ratio as function of relative humidity and temperature*)
Y[phi_?NumericQ, T_?NumericQ] := (rn*phi*psat[T])/(patm - psat[T])

phiPlot = Table[Plot[Y[phi, T], {T, 0, 50}, PlotStyle -> {Thin, Blue}], {phi,0.1, 1, 0.1}];

phiSatPlot = Plot[Y[1, T], {T, 0, 50}, PlotStyle -> Green, PlotRange -> {{0, 50}, {0, 0.050}}];

Tref = 0;

(*YB is humidity ratio as function of enthalpy and temperature*)
YB[h_?NumericQ,T_?NumericQ] := (h - (cpda*(T + 273.15 - (Tref + 273.15))))/((cpwv*(T + 273.15 - (Tref + 273.15))) + hvref)

h1Plot = Table[Plot[YB[h, T], {T, 0, 55}, PlotStyle -> {Thin, Red},PlotRange -> {{0, 50}, {0, 0.050}}], {h, 0, 200, 5}];

Show[h1Plot, phiSatPlot, phiPlot, Frame -> True,
FrameTicks -> {{None, Range[0, 0.050, 0.002]}, {Range[0, 50, 10],
None}}]


So the following graphic shows humidity ratio (Y axis), dry bulb temperature (x axis), relative humidity (blue lines) and enthaply

The issue is that I can't erase the enthalpy (red) lines above the saturation line (green curve). Does anybody knows how I could do that?

• Welcome to the Mathematica Stack Exchange. Could you please include definition for psat?
– Syed
Commented Nov 13, 2022 at 5:12
• phiSatPlot = Plot[Y[1, T], ... What is Y? Please add all this info to your post and delete your comments. Make sure you copy from this web page to a new notebook on a fresh kernel. If you can't see the plot, we can't see it either.
– Syed
Commented Nov 13, 2022 at 10:50
• Thank you @Syed! I did the alterations and tried the code in a new notebook to check the plot. Commented Nov 13, 2022 at 11:29

Clear["Global*"]

c8 = -5800.2206;
c9 = 1.3914993;
c10 = -0.04860239;
c11 = 0.41764768*10^-4;
c12 = -0.14452093*10^-7;
c13 = 6.5459673;
patm = 101325;
Tref = 0;

(*Pa*)(*atm pressure*)rn = 0.622;(*kg/kg*);(*molar ratio*)cp \
= 1.004;(*KJ/kg K*);(*air specific heat*)hvref = \
2501;(*KJ/kg*);(*vapour enthalpy*)cpda = 1.006;(*KJ/kg K*);(*specific \
heat dry air*)cpwv = 1.84;(*KJ/kg K*);(*specific heat water \
vapour*)(*Y is humidity ratio as function of relative humidity and \
temperature*)

Y[phi_?NumericQ, T_?NumericQ] := (rn*phi*psat[T])/(patm - psat[T])

(*psat is saturation pressure calculated by Hyland Wexler correlation*)

psat[T_?NumericQ] :=
Exp[(c8/(T + 273.15)) +
c9 + (c10*(T +
273.15)) + (c11*((T + 273.15)^2)) + (c12*((T +
273.15)^3)) + (c13*Log[T + 273.15])]

(*YB is humidity ratio as function of enthalpy and temperature*)
YB[h_?NumericQ,
T_?NumericQ] := (h - (cpda*(T +
273.15 - (Tref + 273.15))))/((cpwv*(T +
273.15 - (Tref + 273.15))) + hvref)

h1Plot = ListLinePlot[
Table[YB[h, T], {h, 0, 200, 5}, {T, 0, 55}]
, PlotStyle -> {{Thin, Red}}
];

phiPlot =
Plot[Evaluate@Table[Y[phi, T], {phi, 0.1, 0.9, 0.1}], {T, 0, 50}
, PlotStyle -> Blue
, PlotRange -> {0, 0.05}
, PlotRangePadding -> None
, Frame -> True
, FrameTicks -> {
{None, Range[0, 0.050, 0.002]}
, {Range[0, 50, 10], None}
}
, ImageSize -> 500
];

phiSatPlot = Plot[Y[1, T]
, {T, 0, 50}
, PlotStyle -> {Thick, Darker@Green}
, PlotRange -> {{0, 50}, {0, 0.050}}
, FillingStyle -> Directive[Red, HatchFilling[-10 Degree, 0.2, 7]]
, Filling -> {1 -> {Top, Directive[Opacity[1], White]}}

];

Show[phiPlot, h1Plot, phiSatPlot]


• Thank you @Syed! It helped me a lot! Commented Nov 13, 2022 at 19:49

Add a white domain erasePlot to the plot.

erasePlot =
Plot[Y[1, T], {T, 0, 50}, PlotStyle -> Green,
PlotRange -> {{0, 50}, {0, 0.050}}, Filling -> Top,
FillingStyle -> White];
Show[h1Plot, phiSatPlot, phiPlot, erasePlot, Frame -> True,
FrameTicks -> {{None, Range[0, 0.050, 0.002]}, {Range[0, 50, 10],
None}}]


(*psat is saturation pressure calculated by Hyland Wexler correlation*)
psat[T_?NumericQ] :=
Exp[(c8/(T + 273.15)) +
c9 + (c10*(T +
273.15)) + (c11*((T + 273.15)^2)) + (c12*((T +
273.15)^3)) + (c13*Log[T + 273.15])]
c8 = -5800.2206;
c9 = 1.3914993;
c10 = -0.04860239;
c11 = 0.41764768*10^-4;
c12 = -0.14452093*10^-7;
c13 = 6.5459673;
patm = 101325;(*Pa*)(*atm pressure*)
rn = 0.622;(*kg/kg*);(*molar ratio*)
cp = 1.004;(*KJ/kg K*);(*air specific heat*)
hvref = 2501;(*KJ/kg*);(*vapour enthalpy*)
cpda = 1.006;(*KJ/kg K*);(*specific heat dry air*)
cpwv = 1.84;(*KJ/kg K*);(*specific heat water vapour*)

(*Y is humidity ratio as function of relative humidity and \
temperature*)
Y[phi_?NumericQ, T_?NumericQ] := (rn*phi*psat[T])/(patm - psat[T])

phiPlot =
Table[Plot[Y[phi, T], {T, 0, 50}, PlotStyle -> {Thin, Blue}], {phi,
0.1, 1, 0.1}];

phiSatPlot =
Plot[Y[1, T], {T, 0, 50}, PlotStyle -> Green,
PlotRange -> {{0, 50}, {0, 0.050}}];

Tref = 0;

(*YB is humidity ratio as function of enthalpy and temperature*)
YB[h_?NumericQ,
T_?NumericQ] := (h - (cpda*(T +
273.15 - (Tref + 273.15))))/((cpwv*(T +
273.15 - (Tref + 273.15))) + hvref)

h1Plot =
Table[Plot[YB[h, T], {T, 0, 55}, PlotStyle -> {Thin, Red},
PlotRange -> {{0, 50}, {0, 0.050}}], {h, 0, 200, 5}];

erasePlot =
Plot[Y[1, T], {T, 0, 50}, PlotStyle -> Green,
PlotRange -> {{0, 50}, {0, 0.050}}, Filling -> Top,
FillingStyle -> White];
Show[h1Plot, phiSatPlot, phiPlot, erasePlot, Frame -> True,
FrameTicks -> {{None, Range[0, 0.050, 0.002]}, {Range[0, 50, 10],
None}}]

• Thanks @cvgmt! I don´t think I can remove h1Plot entirely since it shows enthalpy lines that are important to the chart. Commented Nov 13, 2022 at 11:57
• @JulioAraujoDosSantos updated. Commented Nov 13, 2022 at 12:07
• thank you very much my friend! It solves the problem. Commented Nov 13, 2022 at 12:56

You can use RegionFunction option with your boundary function Y[1,T]

Note that you don't need to plot each contour separately. The method used below for newh1Plot will save memory.

With definitions as in OP then

newh1Plot =
Plot[
Table[YB[h, T], {h, 0, 200, 5}]
, {T, 0, 55}
, PlotStyle -> {Thin, Red}
, PlotRange -> {{0, 50}, {0, 0.050}}
, RegionFunction -> Function[{x, y}, y <= Y[1, x]]
]


and

Show[newh1Plot, phiSatPlot, phiPlot
, Frame -> True
, FrameTicks -> {{None, Range[0, 0.050, 0.002]}, {Automatic, None}}
]
`

Hope this helps.

• thank you very much, it helps a lot! And also thank you for the memory optmization tip. Commented Nov 13, 2022 at 12:59