Notebook 3.3: Calculating Physical Properties from Thermodynamic Potentials

This notebook shows how physical properties like the isothermal compressibility can be calculated from accurate thermodynamic potentials (in this case, Fiestel’s Gibbs free energy function for ice, used in oceanographic modeling). It also introduces function definition in Maxima.

• Download and install the Euler Math Toolbox.
• Cut and paste the code for notebook3.3 below into a plain text file, and save it as “notebook3.3.en”.
• Cut and paste the code for Gice into a plain text file, and save it as “Gice.e”. This file must be placed in the same directory as notebook3.3.en.
• Double-click notebook3.3.en. It should open up in the Toolbox.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition. Italics mark the items students had to fill in themselves.

--------------------------------------------------notebook 3.3 code begins------------------------------------------
 Notebook 3.3 Calculating Physical Properties from Thermodynamic Potentials
% 
% In this notebook, we'll use thermodynamic potentials to predict material properties.
% In particular, we'll use an empirical Gibbs free energy function for ice to show how 
% properties like density, energy, and entropy vary with temperature and pressure.
% 
% REFERENCES
%     1. Fiestel, R., Wagner, W., Tchijov, V., Guder, C., 
%         "Numerical implementation and oceanographic application of the Gibbs 
%          potential of ice",  Ocean Science, 1, 29-38, 2005.
%     2. Fiestel, R., Wright, D. G., Miyagawa, K., Harvey, A. H, Hruby, A. H., 
%        Jackett, D. R., McDougall, T. J., Wagner, W., "Mutually consistent 
%        thermodynamic potentials for fluid water, ice, and seawater: a new 
%        standard for oceanography", Ocean Science, 4, 275-291, 2008.
% 
% BACKGROUND
%     We can calculate many material properties as first and second derivatives 
%     of the Gibbs free energy (see Section 3.9 in Atkins 8e).
%     We've already seen how functions can be defined in Euler. In this notebook,
%     we'll define functions in Maxima  using the ":=" operator. For example, to 
%     define a function that calculates the area of a circle from its radius,
%     we can type
% 
>: area(r) := %pi*r^2
 
                                                2
                                area(r) := %pi r
 
% 
% The area of a circle with a radius of 2 m can be calculated as
>: area(2)
 
                                      4 %pi
 
% or
>: area(r), r=2
 
                                      4 %pi
 
% 
% ...both mean the same thing, and give the same result.
% 
% We've already seen that we can calculate derivatives like dA/dr 
% with the diff() function:
% 
>: diff(area(r),r)
 
                                     2 %pi r
 
% 
% We can calculate second derivatives as well by repeating the differentiation
% 
>: diff(diff(area(r),r), r)
 
                                      2 %pi
 
% 
% but there's an easier way. Maxima lets you specify the order of the derivative
% as a third argument to diff. 
% For the second derivative of area with respect to r, we can type
% 
>: diff(area(r),r,2)
 
                                      2 %pi
 
% --------------------------------------------------------------------------------------------------------------------------------------------------------------
% In the space below, calculate the surface area of a sphere from its volume. 
% The volume of a sphere is (4/3) pi r^3,  and the surface area of the sphere is 
% dV/dr.
% 
>: V(r) := (4/3)*%pi*r^3
 
                                        4      3
                                V(r) := - %pi r
                                        3
 
>: surfacearea : diff(V(r),r)
 
                                           2
                                    4 %pi r
 
% 
% The line below loads the thermodynamic potential from Reference
% 1 into Maxima from an external file ("Gice.e"). The Maxima function
% Gice(T,P) computes  the Gibbs free energy of one kilogram of ice
% (in J/kg) at a temperature of T (in K) and a pressure of P (in
% Pa).
% 
>load("Gice")
 
% 
% We can have Maxima print back the Gice(T,P) function for us:
>: Gice(T,P)
 
                                                            2
        273.16 (- 65.22370501477501 (1.7162371555667344E-4 T
  + atan2(0.051046477118412, 0.0036608581051398 T + 0.037153909034639)
  (0.0036608581051398 T + 0.037153909034639)
  + atan2(0.051046477118412, 0.037153909034639 - 0.0036608581051398 T)
  (0.037153909034639 - 0.0036608581051398 T)
                                                                    2
  + 0.025523238559206 log((0.0036608581051398 T + 0.037153909034639)
  + 0.0026057428262006) + 0.025523238559206
                                                2
  log((0.037153909034639 - 0.0036608581051398 T)  + 0.0026057428262006)
  + 0.21205793408363) - (- 2.7329787774916599E-11
                                         2
  (0.0016349032219038 P - 165.6565689594)
  + 5.0905901194652602E-5 (0.0016349032219038 P - 165.6565689594)
                                               2
  - 80.98785064626451) (1.9417364295872538E-5 T
  + atan2(0.34331589201784, 0.0036608581051398 T + 0.34509582956282)
  (0.0036608581051398 T + 0.34509582956282)
  + atan2(0.34331589201784, 0.34509582956282 - 0.0036608581051398 T)
  (0.34509582956282 - 0.0036608581051398 T)
                                                                  2
  + 0.17165794600892 log((0.0036608581051398 T + 0.34509582956282)
  + 0.11786580171201) + 0.17165794600892
                                               2
  log((0.34509582956282 - 0.0036608581051398 T)  + 0.11786580171201)
  - 0.045958113581803) + (2.3961751351811602E-11
                                         2
  (0.0016349032219038 P - 165.6565689594)
  - 5.7552976563435302E-5 (0.0016349032219038 P - 165.6565689594)
                                               2
  - 75.8695106343435) (- 1.951803454312407E-5 T
                                                  2
  + (log((0.0036608581051398 T + 0.34509582956282)  + 0.11786580171201)
  (0.0036608581051398 T + 0.34509582956282))/2
                                                  2
  + (log((0.34509582956282 - 0.0036608581051398 T)  + 0.11786580171201)
  (0.34509582956282 - 0.0036608581051398 T))/2
  - 0.34331589201784 atan2(0.34331589201784, 
 0.0036608581051398 T + 0.34509582956282)
  - 0.34331589201784 atan2(0.34331589201784, 
 0.34509582956282 - 0.0036608581051398 T) + 1.034399513744271)
                                              2
  + 45.951447199735 (- 1.2491541583149561E-4 T
                                                   2
  + (log((0.0036608581051398 T + 0.037153909034639)  + 0.0026057428262006)
  (0.0036608581051398 T + 0.037153909034639))/2
                                                   2
  + (log((0.037153909034639 - 0.0036608581051398 T)  + 0.0026057428262006)
  (0.037153909034639 - 0.0036608581051398 T))/2
  - 0.051046477118412 atan2(0.051046477118412, 
 0.0036608581051398 T + 0.037153909034639)
  - 0.051046477118412 atan2(0.051046477118412, 
 0.037153909034639 - 0.0036608581051398 T) + 0.30140605416073))
  + 3333.18160308627 T - 5.7859365867952202E-22
                                         4
  (0.0016349032219038 P - 165.6565689594)
                                                                3
  + 3.40692612753936E-15 (0.0016349032219038 P - 165.6565689594)
                                                               2
  - 1.89952376891314E-8 (0.0016349032219038 P - 165.6565689594)
  + 0.65502999780479 (0.0016349032219038 P - 165.6565689594) - 632578.704355102
 
% 
% We can use this function to compute (and graph) almost all of
% the major state functions for ice. For example, (dG/dP)_T = V,
% so the volume of one kilogram of ice can be calculated using the
% line below.
% 
% Notice that we calculate (dG/dP)_T symbolically with diff(), and
% then evaluate it at t,p using at(). In English, the following
% line says, "V(t,p) is the value of dG/dP at temperature t and
% pressure p".
% 
% 
>: V(t,p) := at(diff(Gice(T,P), P), [T=t, P=p])
 
               V(t, p) := at(diff(Gice(T, P), P), [T = t, P = p])
 
% 
% At 263.15 K and 10101325 Pascals, the volume of one kilogram of ice (in 
% cubic meters) is
>: V(263.15,10101325)
 
                               0.0010878711964292
 
% 
% Let's plot the volume of 1 kilogram of ice as a function of temperature
% from 250 K to 280 K, at a constant pressure of 1 atm. 
% 
% First, we define lists of temperatures and volumes to plot. It's
% going to be easier if the temperature list can be read by both
% Maxima and Euler. Notice that
% 
% >Temperature = 250:280  would define the array ONLY in Euler,
%%>Temperature ::= 250:280 defines the array in BOTH Maxima and
% Euler. 
% 
% If the ::= operator doesn't work, you're using an old version
% of Euler. Upgrade.
> Temperature ::=250:280; 
% 
% Now we use the Maxima function V(t,p) to make a list of volumes
% in Euler. Recall that the mxmeval function sends a command (in
% double quotes) from Euler to Maxima.
> Volume = mxmeval("V(Temperature,p)", p=101325);
> plot2d(Temperature,Volume);
> xlabel("Temperature (K)");
> ylabel("Volume (m^3/kg)");
> title("Volume of 1 kg ice, predicted by Fiestel's 2005 Gibbs Potential");
> insimg;

>Volume = mxmeval("V(Temperature,p)", p=101325)*(1/1e-6)*(1/1000);
>plot2d(Temperature, Volume);
>label("1 atm",266,mean(Volume));
>Volume = mxmeval("V(Temperature,p)", p=500*101325)*(1/1e-6)*(1/1000);
>plot2d(Temperature, Volume, add=1, thickness=2);
>label("500 atm",266,mean(Volume));
>Volume = mxmeval("V(Temperature,p)", p=1000*101325)*(1/1e-6)*(1/1000);
>plot2d(Temperature, Volume, add=1, thickness=2);
>label("1000 atm",266,mean(Volume));
> xlabel("Temperature (K)");
> ylabel("Volume (cm^3/g)");
> title("Specific Volume of Ice Predicted by Fiestel's 2005 Gibbs Potential");
>insimg;

% ----------------------------------------------------------
%%In the space below, define a function that calculates the coefficient
% of isothermal compressibility
% 
% kappa_T = -(1/V)(dV/dP)_T
% 
% for ice as a function of temperature and pressure. Call your function
% to calculate the isothermal compressibility of ice at the triple
% point; under those conditions, kappa_T is about 1.18 x 10^-10
% Pa^-1.
% 
% 
>: kappaT(t,p) := at(-(1/V(T,P))*diff(V(T,P),P),[T=t,P=p]);
 
                                 1
        kappaT(t, p) := at((- -------) diff(V(T, P), P), [T = t, P = p])
                              V(T, P)
 
>: kappaT(273.16,611.657)
 
                             1.1791917359556073E-10
 
% --------------------------------------------------------------
%%In the space below, plot the isothermal compressibility for ice
% predicted for temperatures from 250 K to 280 K, at pressures of
%  1, 100, and 500 atm. Scale the y-axis of the plot so it is easy
% to read, and label each individual curve on your plot.
% 
>scalefactor=10^floor(log(mxmeval("kappaT(t,p)",t=mean(Temperature),p=101325))/log(10))
              1e-011 
>kappaTvalues = mxmeval("kappaT(t,p)",t=Temperature,p=101325)/scalefactor;
>plot2d(Temperature,kappaTvalues);
>xlabel("Temperature (K)");
>ylabel("Isothermal Compressibility kappa_T (Pa^-1)/" | scalefactor);
>title("Isothermal Compressibility of Ice Predicted by Fiestel's 2005 Gibbs Potential");
>label("1 atm",mean(Temperature),mean(kappaTvalues));
>kappaTvalues = mxmeval("kappaT(t,p)",t=Temperature,p=100*101325)/scalefactor;
>plot2d(Temperature,kappaTvalues,add=1);
>label("100 atm",mean(Temperature),mean(kappaTvalues));
>kappaTvalues = mxmeval("kappaT(t,p)",t=Temperature,p=500*101325)/scalefactor;
>plot2d(Temperature,kappaTvalues,add=1);
>label("500 atm",255,0.85);
>insimg;

-------------------------------notebook 3.3 code ends------------------------------------------------------------------

-------------------------------Gice code begins------------------------------------------------------------------------------
:Gice(T,P) :=  block([g0: [-632578.704355102, 0.655029997804786, -1.89952376891314e-8,
 3.40692612753936e-15,-5.78593658679522e-22],	t: [3.71539090346389e-2 + 
5.10464771184122e-2*%i, 0.345095829562823 + 0.343315892017841*%i], r: 
[45.951447199735 + 65.223705014775*%i, 0.0 + 0.0*%i],r2: [-75.8695106343435 
-80.9878506462645*%i, -5.75529765634353e-5 + 5.09059011946526e-5*%i, 
2.39617513518116e-11 -2.73297877749166e-11*%i], Tt: 273.160, Pt: 611.657, p0: 
165.6565689594, Pr: P/Pt, Tr: T/Tt,  sigma0: -3333.18160308627], r[2]: sum(r2[k+1]*(P/Pt-
p0)^k,k,0,2), -sigma0*T + sum(g0[k+1]*(P/Pt-p0)^k,k,0,4)+Tt*realpart(sum(r[k]*((t[k]-T
/Tt)*log(t[k]-T/Tt) + (t[k]+T/Tt)*log(t[k]+T/Tt) - 2*t[k]*log(t[k])-(T/Tt)^2/t[k]),k,1,2)))$
% 
function atan2(y,x) := atan(y/x)
------------------------------Gice code ends--------------------------------------------------------------------------------

Notebook 3.2: Predicting Boiling Points from Liquid/Vapor Gibbs Free Energy Functions

This notebook shows how Shomate functions for a liquid and vapor can be used to predict boiling points. It also introduces function programming and for loops.

* Download and install the Euler Math Toolbox.
* Cut and paste the code below into a plain text file.
* Save the file with an .en file extension.
* Double-click the file. It should open up in the Toolbox.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition. Italics mark the items students had to fill in themselves.

----------------------------------------------------snip snip----------------------------------------------------------
 Notebook 3.2: Variation of the Gibbs Free Energy with Temperature
% 
% In this notebook, we'll predict the boiling points of a liquid by looking for a 
% temperature that makes the Gibbs free energies of the vapor and liquid equal. 
% To compute the Gibbs free energy at standard state (constant pressure), 
% we'll integrate the equation
% 
%               dG = -SdT
% 
% We'll need to find S as a function of T to do the integral. A very accurate 
% empirical fit called the Shomate equation can be used here:
% 
% S = A*log(t) + B*t + C*t^2/2 + D*t^3/3 - E/(2*t^2) + G;
% 
% where we can find look up the Shomate coefficients A, B, C, D, E, and G on 
% the NIST Chemistry Webbook (http://webbook.nist.gov/chemistry). 
% 
% When we're doing a complicated calculation that will be repeated many times, 
% it's convenient to define the calculation as a "function". In Euler, we can define
% a function with a series of statements like
% 
%      function f(x)
%          ....
%      return y
%      endfunction
% 
% where f is the function name, x is value passed to the function, and y is the
% calculated value of f(x) to be returned. For example, here is a function that 
% computes the pressure of an ideal gas from its temperature and molar volume:
%
>function Pideal(V,T)
$R=0.082059;
$return R*T/V;
$endfunction
% 
% We can compute the pressure of an ideal gas with V = 22.4 L/mol and 
% T = 273.15 K in atmospheres as
>Pideal(22.4, 273.15)
      1.000643564732 
% 
% which is exactly the same as typing
>0.082059*273.15/22.4
      1.000643564732 
% 
% Let's define a function that computes the entropy of steam at temperature 
% T using the Shomate equation. The Shomate coefficients for water and steam
% are taken from 
% http://webbook.nist.gov/cgi/cbook.cgi?ID=C7732185&Units=SI&Mask=7.
% 
>function EntropySteam(T) 
$t=T/1000;
$A=30.09200;
$B=6.832514;
$C=6.793435;
$D=-2.534480;
$E=0.082139;
$G=223.3967;
$return A*log(t) + B*t + C*t^2/2 + D*t^3/3 - E/(2*t^2) + G;
$endfunction
% 
% This function computes the entropy in SI units (J mol^-1 K^-1 ) 
% To check to see that we've typed all the coefficients in correctly, we can 
% compare our entropy at 298.15 with the CODATA Review value, 
% 188.835 +/- 0.010 J mol^-1 K^-1 
% (Cox, J. D.; Wagman, D. D.; Medvedev, V. A., CODATA Key Values for 
% Thermodynamics, Hemisphere Publishing Corp., New York, 1984, 1). 
% 
>EntropySteam(298.15)
      188.8352691656 
% 
% Looks good. Now let's calculate the entropy of liquid water, again using 
% Shomate equation coefficients from the NIST Webbook:
% 
>function EntropyWater(T) 
$t = T/1000;
$A = -203.6060;
$B = 1523.290;
$C = -3196.413;
$D = 2474.455;
$E = 3.855326;
$G = -488.7163;
$return A*log(t) + B*t + C*t^2/2 + D*t^3/3 - E/(2*t^2) + G;
$endfunction
% 
% Again, we see that the function correctly predicts the CODATA Review value
%  of the standard entropy of water at
% 298.15 K, which should be 69.95 ± 0.03 J mol^-1 K^-1.
% 
>EntropyWater(298.15)
      69.95364324026 
% 
% To integrate -S dT from 298.15 K to 500 K, we can type
% 
>-integrate("EntropySteam(x)",298.15,500)
     -40041.68164763 
% 
% Typing "EntropySteam(T)" won't work; because Euler doesn't know we want to 
% integrate with respect to T. 
% We have to call the integration variable "x".  
% 
% The integral gives us the change in Gibbs free energy between 298.15 and 
% 500 K. To get the Gibbs free energy itself at some temperature T, we can 
% add the Gibbs free energy of formation at 298.15 to the integral of -S dT
% from 298.15 K to T:
% 
>function GSteam(T)
$DeltaGfo = -228570;
$return DeltaGfo - integrate("EntropySteam(x)",298.15,T)
$endfunction
% 
% At 500 K, the Gibbs free energy of steam is
% 
>GSteam(500)
     -268611.6816476 
% 
% Now let's build an array of Gibbs free energies for the gas (call it Ggas)
% and a corresponding array of temperatures. Suppose we want 200 points:
>Ggas = 1:200; T = 1:200;
% 
% Unfortunately, we can't just type Ggas = GSteam(T), because 
% our GSteam function only calculates one Gibbs free energy at a time. We
% can use a "for" loop to process each individual element of the T and Ggas 
% arrays. The syntax is
% 
%     for index=start to finish;  ... end
% 
% For example, if we wanted to set up 200 temperature points, starting from
% 299.15 K and ending with 499.15 K, we could type
% 
>for i=1 to 200; T[i] = 298.15 + i; end
% 
% Now set up the Ggas array. We'd like the Gibbs free energy to
% be displayed in kilojoules per mole.
% 
>for i=1 to 200; Ggas[i]= GSteam(T[i])/1000; end
% 
% Plot the Gibbs free energy of the gas as a function of T:
% 
>plot2d(T,Ggas,thickness=2);
>xlabel("Temperature (K)");
>ylabel("Gibbs Free Energy (kJ/mol)");
>label("gas",T[160],Ggas[150]);
% 
% Now let's add a liquid curve to the plot. First, define a function
% to calculate the Gibbs free energy of water at a specific temperature:
% 
>function GWater(T)
$DeltaGo = -237130;
$return DeltaGo -integrate("EntropyWater(x)",298.15,T)
$endfunction
% 
% Now set up an array to hold the Gibbs free energy data for the liquid
% (call it Gliq). 
% 
>Gliq = 1:200;
>for i=1 to 200; Gliq[i] = GWater(T[i])/1000; end
% 
% Add the liquid curve to the graph.
% 
>plot2d(T,Gliq,add=1,thickness=1);
>title("Determination of Boiling Point from Vapor/Liquid Gibbs Free Energies");
>label("liquid",T[160],Gliq[150]);
% 
% When the Gibbs free energies of the liquid and gas are equal,
% the liquid boils. To get the boiling point Tbp, then, we need
% to solve the equation GSteam(Tbp) = GWater(Tbp) for Tbp.
% 
% In a previous notebook, we used Euler's root() function to solve
% equations like this. We need an initial guess for Tbp. Anything
% that isn't too far off will do, so we'll use the midpoint of our
% temperature range as a rough guess:
>bp = mean(T)
              398.65 
% 
% Now use root to find the intersection of the two curves. 
>bp = root("GSteam(bp)-GWater(bp)",bp)
       373.216970275 
% 
% ...which is only about .05 K different from the true standard
% boiling point of water.
% 
% We can double-check by printing the Gibbs free energies of the
% liquid and gas at this temperature. They must be equal.
>GSteam(bp)
     -243040.0752017 
>GWater(bp)
     -243040.0752017 
% 
% Finally, let's mark the boiling point on the plot. We want to
% drop a line from the intersection at (bp, GSteam(bp)/1000) to
% the x axis at (bp, min(Ggas)). 
%% 
>plot2d(bp,GSteam(bp)/1000,add=1,points="[]");
>plot2d([bp,bp],[GSteam(bp)/1000,min(Ggas)],add=1,style="--");
>label("Boiling Point",bp,Ggas[190]);
>insimg;

% 
% Make a NEW notebook (using this one as a model) that locates the
% 
% boiling point of some other compound by plotting the Gibbs free
% energy for the phases in equilibrium. Your notebook should include
% a graph set up and labeled like the one above. It should also
% include a numerical estimate of the boiling or melting temperature.
% 
>

How to set styles on text input fields in ActionScript 3.0

Problem: Setting the text style on an input TextField doesn’t work as expected; for example,


...
var t:TextField = new TextField();
var f:TextFormat = new TextFormat();
...
t.type = TextFieldType.INPUT;
f.font = "Arial";
f.size = 20;
t.setTextFormat(f);
addChild(t);
...

Text typed into the text field appears in the default serif font, not in Arial 20.

Solution You can’t style user input or replaced text in Actionscript 3.0 with setTextFormat, instead, you have to use defaultTextFormat. Just replace the line that calls setTextFormat with


t.defaultTextFormat = f;


and the input field will use the text format in f correctly. This won’t work if you’re using an ActionScript 3.0 style sheet with the textfield, though.

Notebook 2.1: Work in Reversible Isothermal Expansions

This short exercise gives students more practice using Maxima’s integration and comma operators.

• Download and install the Euler Math Toolbox.
• Cut and paste the code below into a plain text file.
• Save the file with an .en file extension.
• Double-click the file. It should open up in the Toolbox.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition. Italics mark the items students had to fill in themselves.

--------------------------------------------------------snip here----------------------------------------------------
 Notebook 2.1: Work in Reversible Isothermal Expansions
% 
% In reversible expansions or compressions, the external pressure constantly changes so that 
% it remains exactly equal to the pressure of the gas. In this special case, we can substitute
% the pressure calculated from the gas's equation of state for the external pressure when calculating work.
% For example, the work done to compress n moles of an ideal gas from V1 to V2 under reversible
% conditions at temperature T is
% 
>: assume(V1>0, V2>0, V1>V2, V>0)
 
                        [V1 > 0, V2 > 0, V1 > V2, V > 0]
 
>: work_ideal: -integrate(n*R*T/V,V,V1,V2)
 
                            n R T (log(V1) - log(V2))
 
% If the "assume" had been omitted, Maxima would have asked a series of questions about 
% the signs of the volumes during the integration.
% 
% If we want to calculate the work for specific temperatures and volumes, we can use Maxima's
% comma operator. For example,
>: work_ideal, n=1, R=8.314, T=298, V1=1, V2=0.5, numer
 
                                1717.322046434265
 
% 
% Integrating -n*R*T/V with respect to V from V1 to V2 at constant T gives -n*R*T*ln(V2/V1), so we could
% also have calculated the work as -n*R*T*log(V2/V1). 
%  --> Notice that the Maxima function for natural logarithms is log(),  NOT ln()!)
%  --> The final ",numer" tells Maxima that we want a number as an answer. Maxima will sometimes 
%        just display logs of whole numbers without calculating them in a final result if we don't include this.
% ---------------------------------------------------------------------------------------------------------
% In the space below, create an expression called work_vdw that calculates the work for a reversible isothermal 
% expansion of a van der Waals gas as it expands from volume V1 to volume V2 at temperature T.
% 
>: assume(V>b, V1>b, V2>b)
 
                             [V > b, V1 > b, V2 > b]
 
>: work_vdw: -integrate(n*R*T/(V-b)-n^2*a/V^2,V,V1,V2)
 
                                       2                          2
                                    a n                        a n
              - n R T log(V2 - b) - ---- + n R T log(V1 - b) + ----
                                     V2                         V1
 
% ---------------------------------------------------------------------------------------------------------
% In the space below, use the comma operator to calculate the work for one mole of  xenon expanding from 1 m^3
% to a volume of 10 m^3 at a temperature of 298 K.
% 
>: work_vdw, a=5.125, b=1.06e-2, T=298, n=1, V1=1, V2=10, R=8.314,numer
 
                               - 5723.982679505956
 
% ---------------------------------------------------------------------------------------------------------
% In the space below, use the comma operator to calculate the work for one mole of  xenon expanding from 1 m^3
% to a volume of 10 m^3 at a temperature of 298 K IGNORING INTERMOLECULAR ATTRACTIONS.
% 
>: work_vdw, a=0, b=1.06e-2, T=298, n=1, V1=1, V2=10, R=8.314, numer
 
                               - 5728.595179505956
 
% ---------------------------------------------------------------------------------------------------------
% In the space below, use the comma operator to calculate the work for one mole of  xenon expanding from 1 m^3
% to a volume of 10 m^3 at a temperature of 298 K IGNORING INTERMOLECULAR REPULSIONS.
% 
>: work_vdw, a=5.125, b=0, T=298, n=1, V1=1, V2=10, R=8.314, numer
 
                               - 5700.207854019444
 
% ---------------------------------------------------------------------------------------------------------
% In the space below, use the comma operator to calculate the work for one mole of  xenon expanding from 1 m^3
% to a volume of 10 m^3 at a temperature of 298 K IGNORING ALL INTERMOLECULAR FORCES.
% 
>: work_vdw, a=0, b=0, T=298, n=1, V1=1, V2=10, R=8.314, numer
 
                               - 5704.820354019445

Notebook 2.0: Introducing integration with Maxima

This very simple Euler notebook introduces integration with Maxima by asking students to set up simple integrals for calculating expansion work.

To use the notebook:

• Download and install the Euler Math Toolbox.
• Cut and paste the code below into a plain text file.
• Save the file with an .en file extension.
• Double-click the file. It should open up in the Toolbox.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition. Italics mark the items students had to fill in themselves.

--------------------------------------------------------snip here----------------------------------------------------

 Notebook 2.0 - Integration with Maxima
% 
% When we're finding derivatives, we are finding the instantaneous rate
% of change for a function. For example, we can see how P changes with 
% V for an ideal gas by calculating the derivative dP/dV at constant T:
 
>: Pideal : R*T/V
>: dPdV_ideal : diff(Pideal,V)
 
                                        R T
                                      - ---
                                         2
                                        V
 
>: Pvdw : R*T/(V-b)-a/V^2
>: dPdV_vdw : diff(Pvdw, V)
 
                                 2 a     R T
                                 --- - --------
                                  3           2
                                 V     (V - b)
 
% 
% Often we want to calculate the change in a variable from its rate
% of change. For example, if we have dP/dV, we can reverse the process
% above to find P, by integrating the derivative with respect to
% V:
% 
>: Pideal : integrate(dPdV_ideal,V)
 
                                       R T
                                       ---
                                        V
 
>: Pvdw : integrate(dPdV_vdw,V)
 
                                    R T    a
                                   ----- - --
                                   V - b    2
                                           V
 
% 
% Integration undoes differentiation, and vice versa:
% 
>: integrate(diff(Pideal,V),V)
 
                                       R T
                                       ---
                                        V
 
>: diff(integrate(Pideal, V),V)
 
                                       R T
                                       ---
                                        V
 
% -------------------------------------------------------------------------------------
% In the space below, calculate the derivative of P for a Dieterici gas with respect
% to V (with T held constant. Call it dPdV_Dieterici.
% 
 
>: P_Dieterici : R*T*exp(-a/(R*T*V))/(V-b)
>: dPdV_Dieterici : diff(P_Dieterici,V)
 
% -------------------------------------------------------------------------------------
% In the space below, show that integrating dPdV_Dieterici gives back
% P for a Dieterici gas..
% 
>: integrate(dPdV_Dieterici,V)
% -------------------------------------------------------------------------------------
% We can use integration to calculate the work done during gas expansions, or the work
% required to compress a gas. This is done by integrating the external pressure (NOT the
% gas pressure) with respect to volume, over some volume range. 
% 
% In Maxima, you can integrate over a range by adding the beginning and end of the range as 
% third and fourth arguments for the integrate() function. For example, the work done to compress a gas from 4 m^3 to 2 m^3 under an external
% pressure of 101000 Pa is
 
>: work: -integrate(101000,V,4,2)
 
% When Pext is constant, the work is simply -Pext times the volume change, so we could have
% calculated the work more simply as
% 
>: work: -101000*(2-4)
% 
% Both calculations give work in joules, because we've used SI units (pascals and cubic meters)
% throughout the calculation. 
% 
% In the space below, calculate the work done (in joules) when a sample of 1.70 mol of argon gas
% expands from 20 L to 50 L against a constant external pressure equal to the final
% pressure of the gas. The temperature of the gas is 25 degrees C.
% 
>: n: 1.70
>: R: 8.314
>: T: 25+273.15
>: V2: 50*1e-3
>: V1: 20*1e-3
>: Pext: n*R*T/V2
>: work: -integrate(Pext,V,V1,V2)
% or, more simply,
>: work: -Pext*(V2-V1)
% -------------------------------------------------------------------------------------
% In the space below, calculate the expansion work done when 1 mole of water is electrolysed
% under a constant 1 atm external pressure at 298 K. Hint: Study example 2.1 in the text!
% 
% Solution: The work is equal to -Pext times the volume change. The initial volume of the water
% is 1 mol * 18 g/mol * 1 mL / 1 g * 1e-3 L/ 1 mL * 1 dm^3 / 1 L * 1e-3 m^3 / 1 dm^3
>:  V1: 1*18*1*1e-3*1*1e-3
% 1/2 mol of O2 and 1 mole of H2 is produced by the reaction. The final pressure of the gases
% will be equal to the external pressure. If the gases are assumed to be ideal, the final volume is
% n*R*T/Pext.
>: n: 1.5
>: R: 8.314
>: T: 298
>: Pext: 1 * 101325
>: V2: n*R*T/Pext
% and the work is
>: work: -integrate(Pext,V,V1,V2)
% or, equivalently, 
>: work: -Pext*(V2-V1)
% We could have simplified the solution by neglecting V1, since gases have about 1000 times
% the molar volume of liquids around room temperature and pressure. Then the work is 
% approximately
>: work: -Pext*V2
------------------------------------snip here------------------------------------

Notebook 1.5: Critical Point Coordinates from Equations of State

This notebook uses Euler’s symbolic math capabilities to derive the critical point pressure, volume, temperature, and compression factor for a van der Waals gas. It then asks students to repeat the calculations for a Dieterici gas.

To use the notebook:

• Download and install the Euler Math Toolbox.
• Cut and paste the code below into a plain text file.
• Save the file with an .en file extension.
• Double-click the file. It should open up in the Toolbox.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition. Italics mark the items students had to fill in themselves.

------------------------------------snip here-----------------------
 Notebook 1.5: Calculating Critical Point Coordinates from Equations of State
% 
% In this notebook, we'll use Euler's symbolic math capabilities to
% derive the critical point pressure, volume,
% temperature, and compression factor for a van der Waals gas. You'll
% then repeat the calculations yourself for a Dieterici gas.
% 
% We saw in the last notebook that Euler calls a separate package 
% called Maxima to do symbolic math.
% When you want numerical answers, use Euler; but if you want to 
% solve or manipulate equations, send your commands to
% Maxima by typing a colon (':') after the prompt.
% 
% Here is the van der Waals equation, in Maxima:
 
>: P : R*T/(V-b)-a/V^2
 
                                    R T    a
                                   ----- - --
                                   V - b    2
                                           V
 
% 
% Notice that Maxima uses a colon rather than an equals sign to set
% variables. Notice also that the result is symbolic, rather than 
% numerical, so you don't have to assign numbers to R, T, a, etc. as 
% you must from a regular Euler prompt.
% 
% Maxima doesn't simplify expressions automatically. You have to be
% explicit about how you want an expression 'simplified'. For example,
% you can remove parent
 
% Maxima can differentiate expressions with the diff(expr,x) function, 
% where expr is the expression and x is the variable you want to 
% differentiate. Let's calculate the first derivative of P with respect
% to V, and call it dPdV:
% 
>: dPdV : diff(P,V)
 
                                 2 a     R T
                                 --- - --------
                                  3           2
                                 V     (V - b)
 
% 
% Call the second derivative d2PdV2:
% 
>: d2PdV2 : diff(dPdV,V)
 
                                  2 R T     6 a
                                 -------- - ---
                                        3    4
                                 (V - b)    V
 
% 
% We could have calculated the second derivative more simply as
% 
>: diff(P,V,2)
 
                                  2 R T     6 a
                                 -------- - ---
                                        3    4
                                 (V - b)    V
 
% 
% where the third argument to diff is the order of the derivative.
% 
% Now let's calculate the temperature and volume of a van der Waals
% gas at the critical point, using Maxima's solve function.
% The first argument of solve is an array of equations to solve.
% The second argument is an array of unknowns to solve for.
% 
% 
>: solve([dPdV=0,d2PdV2=0],[T,V])
 
                                            8 a
                     [[T = 0, V = b], [T = ------, V = 3 b]]
                                           27 b R
 
% 
% We could also have typed 
% 
>: solve([diff(P,V)=0, diff(P,V,2)=0], [T, V])
% 
% for the same result. We can't put P into the array of unknowns,
% because we've assigned an expression to P already.
% 
% Notice that there are two solutions- one at absolute zero, where all
% of the volume is taken up by the molecules- and the "real" critical 
% point. Let's call the second solution "criticalpoint":
>: criticalpoint : %[2]
 
                                    8 a
                              [T = ------, V = 3 b]
                                   27 b R
 
% 
% The % alone means "the last thing calculated in Maxima". The last
% thing was an array, and we are selecting the second element with the 
% bracketed 2.
% 
% We can now use Maxima's "at" function to calculate the coordinates
% of the critical point. The at function's first argument is the 
% variable we want to solve for; the second argument is an array of
% equations that hold at the point where the solution exists.
% 
>: Pc : at(P,criticalpoint)
 
                                        a
                                      -----
                                          2
                                      27 b
 
>: Tc : at(T,criticalpoint)
 
                                      8 a
                                     ------
                                     27 b R
 
>: Vc : at(V,criticalpoint)
 
                                       3 b
 
% 
% Finally, we calculate Zc, the compressibility factor at the critical
% point:
% 
>: Zc : Pc*Vc/(R*Tc)
 
                                        3
                                        -
                                        8
 
% 
% -----------------------------------------------------------------
% In the space below, derive expressions for Pc, Vc, Tc, and Zc for
% a Dieterici gas. You can check your results against Table 1.7 in
% Atkins. 
% 
% We'll clear all the variables above using 'reset':
>reset;
>: P : R*T*exp(-a/(R*T*V))/(V-b)
 
                                            a
                                        - -----
                                          R T V
                                  R T %e
                                  -------------
                                      V - b
 
>: dPdV : diff(P,V)
 
                                   a               a
                               - -----         - -----
                                 R T V           R T V
                           a %e          R T %e
                           ----------- - -------------
                            2                     2
                           V  (V - b)      (V - b)
 
>: d2PdV2 : diff(dPdV, V)
 
                    a               a                a                 a
                - -----         - -----          - -----           - -----
                  R T V     2     R T V            R T V             R T V
          2 a %e           a  %e           2 a %e          2 R T %e
        - ------------- + -------------- - ------------- + ---------------
            3                  4             2        2              3
           V  (V - b)     R T V  (V - b)    V  (V - b)        (V - b)
 
>: solve([dPdV=0,d2PdV2=0],[T,V])
 
                                             a
                     [[T = 0, V = b], [T = -----, V = 2 b]]
                                           4 b R
 
>: criticalpoint : %[2]
 
                                     a
                              [T = -----, V = 2 b]
                                   4 b R
 
>: Pc : at(P,criticalpoint)
 
                                       - 2
                                     %e    a
                                     -------
                                         2
                                      4 b
 
>: Tc : at(T,criticalpoint)
 
                                        a
                                      -----
                                      4 b R
 
>: Vc : at(V,criticalpoint)
 
                                       2 b
 
>: Zc : Pc*Vc/(R*Tc)
 
                                         - 2
                                     2 %e
 
% 
% Is Zc for the Dieterici equation different than Zc for the van der
% Waals eqution? Which is closer to the experimental values of Zc given
% in your text?
% 
>3/8
                   0.375 
>2*%e^-2
         0.2706705664732 
% ANSWER: The experimental values in Table 1.5 range from 0.27 to about 
% 0.31; The Dieterici equation clearly gives more realistic results than
% the van der Waals equation.
------------------------------------snip here----------------

Notebook 1.4: Trends in Boyle temperatures

All of those ugly, clunky-looking equations of state for real gases have to simplify to PV = nRT at very low pressure. It’s tempting to say that real gases behave like ideal gases under those conditions, but that isn’t quite right.

Even at pressures near zero, all of the properties of real gases won’t become ideal. In particular, any property that depends on the second (and higher) virial coefficient won’t necessarily converge to its ideal gas value as pressure drops. That convergence occurs only at the Boyle temperature, where the second virial coefficient is zero …as long as we assume all the higher order virial coefficients are zero, too. Usually that isn’t a bad assumption.

The Euler notebook posted below asks students to find the Boyle temperatures of noble gases from the van der Waals equation of state. They then critically examine trends in their calculated Boyle temperatures. The notebook introduces algebraic calculations with Maxima, including limits, differentiation, and equation solving.

To use the notebook:

• Download and install the Euler Math Toolbox.
• Cut and paste the code below into a plain text file.
• Save the file with an .en file extension.
• Double-click the file. It should open up in the Toolbox.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition. Italics mark the items students had to fill in themselves.

—————————-snip here———————————————
 Notebook 1.4: Boyle Temperatures of Noble Gases
% 
% In this notebook, we'll use Euler's symbolic math capabilities to
% derive an expression for the Boyle temperature in terms of the 
% van der Waals constants a and b. You'll then use this expression
% to try to discover periodic trends in the Boyle temperature.
% 
% We can derive algebraic expressions in a notebook using Euler.
% Symbolic algebra and calculus can be done with a software package 
% called Maxima which is built into Euler. We can talk to Maxima
% directly by typing a colon (':') after the '>' prompt.
% For example, here is the van der Waals equation in Maxima:
% 
>: (P+a/V^2)*(V-b)=R*T
 
                              a
                             (-- + P) (V - b) = R T
                               2
                              V
 
% 
% If you forget to type : after the > prompt, the command gets
% sent to Euler, not to Maxima. Euler will give you an error 
% message (because it can't compute a result unless you assign numbers
% to each variable). 
% 
% It's useful to name equations, variables, and expressions so
% we don't have to type them over and over again. Maxima lets us
% type a name (say, 'vdw') followed by a colon to name an 
% expression: 
% 
>: vdw : (P+a/V^2)*(V-b)=R*T
 
                              a
                             (-- + P) (V - b) = R T
                               2
                              V
 
% 
% Use a colon to set a name, not an equal sign- you can see 
% from the previous line that '=' means equality, not assignment!
% 
% Now we can just type vdw when we want the van der Waals equation:
% 
>: vdw
 
                              a
                             (-- + P) (V - b) = R T
                               2
                              V
 
% 
% We can make substitutions in the equation using the at() function,
% which looks like at(expression,[conditions]). For example, to see
% that the van der Waals equation is the ideal gas 
% equation when a and b are zero, we could type
% 
>: at(vdw,[a=0,b=0])
 
                                    P V = R T
 
% 
% To solve an equation in Maxima, type solve(equation,x), where
% equation is the equation you want to solve and x is the variable
% to solve for. For example, to solve the equation 
% 2*x^2 + 3*x + 1 = 0 for x, you could type 
% 
>: solve(2*x^2 + 3*x + 1 = 0, x)
 
                                      1
                               [x = - -, x = - 1]
                                      2
 
% 
% Notice that the result is an array; there are two 
% solutions to this quadratic equation.  
% 
% To solve the van der Waals equation for P, you could type
% 
>: solve(vdw,P)
 
                                      2
                                 R T V  - a V + a b
                            [P = ------------------]
                                      3      2
                                     V  - b V
 
% 
% Notice that solve returns an array of solutions - even when 
% there is only one solution, as is the case here. We can 
% extract the % solution from the array by putting a subscript
% 1 (which is written as [1]) after solve:
% 
>: solve(vdw,P)[1]
 
                                      2
                                 R T V  - a V + a b
                             P = ------------------
                                      3      2
                                     V  - b V
 
% 
% The square brackets are gone. Suppose we're just interested
% in the right hand side of this equation. We can use 
% Maxima's rhs function to extract the right hand side:
% 
>: Pvdw : rhs(solve(vdw,P)[1])
 
                                    2
                               R T V  - a V + a b
                               ------------------
                                    3      2
                                   V  - b V
 
% 
% This doesn't look like equation 1.21a from the textbook-
% but it is! You can check to see that Pvdw really is 
% equivalent to R*T/(V-b) - a/V^2 by typing
% 
>: is(equal(Pvdw, R*T/(V-b) - a/V^2))
 
                                      true
 
% 
% Maxima reports and uses expressions in the form of rational polynomials
% (ratios of polynomials). As physical chemists, we would rather
% write expressions broken into separate terms that have some physical
% meaning.
% 
% Maxima can rearrange or simplify the expression, but we have to
% tell it what we want. Here, the partfrac() function can be used.
% The call looks like partfrac(expression,x) where the expression
% is expanded into fractions with powers of x in the denominator.
% Here we can use
% 
>: partfrac(Pvdw,V)
 
                                    R T    a
                                   ----- - --
                                   V - b    2
                                           V
 
% 
% which is exactly how the textbook writes the van der Waals
% equation.
% 
% We can also check to see if Pvdw is correct by seeing that Pvdw
% becomes equal to Pideal as V approaches infinity (which is typed
% as 'inf' in Maxima). To calculate the limit of y as x approaches
% a, type limit(y,x,a). If you enter
>: limit(Pvdw,V,inf)
 
 Is  R T  positive, negative, or zero?
 
>: positive
 
                                        0
 
% Maxima asks whether RT is positive, negative, or zero, and
% you must answer to get the limit. (When you answer, type : after
% the prompt, or the answer goes to Euler, not to Maxima) 
% Maxima then tells us that 
% Pvdw approaches 0 as V approaches infinity, which isn't the same
% as saying it approaches Pideal. Instead, we can define Pideal
% and see whether Pideal/Pvdw approaches 1 as V approaches infinity:
% 
>: Pideal : R*T/V
 
                                       R T
                                       ---
                                        V
 
>: limit(Pideal/Pvdw,V,inf)
 
                                        1
 
% We can also check derivatives of Pvdw this way. As Pvdw approaches
% Pideal, all of the derivatives of Pvdw should approach corresponding
% derivatives of Pideal. 
% 
% It's easy to calculate derivatives in Maxima. The derivative of
%%y with respect to x would be typed diff(y,x). To calculate the
% derivative of Pvdw with respect to V, we can type
>: dPdV_vdw : diff(Pvdw,V)
 
                                  2                2
                2 R T V - a   (3 V  - 2 b V) (R T V  - a V + a b)
                ----------- - -----------------------------------
                  3      2                 3      2 2
                 V  - b V                (V  - b V )
 
% In the space below, calculate the derivative of Pideal with respect
% to V and name it dPdV_ideal.
% 
>: dPdV_ideal : diff(Pideal,V)
 
                                        R T
                                      - ---
                                         2
                                        V
 
% In the space below, show that dP/dV for a van der Waals gas approaches
% dP/dV for an ideal gas as V approaches infinity.
% 
>: limit(dPdV_ideal/dPdV_vdw,V,inf)
 
                                        1
 
% 
% ----------------------------------------------------------------
% 
% Now let's derive an exptression for the Boyle temperature of a
% van der Waals gas. Our approach will be to:
%%1. Calculate dZ/d(1/V) for a van der Waals gas.
%%2. Calculate B as the limit of dZ/d(1/V) as V approaches
%%infinity.
%%3. Find the Boyle temperature as the the temperature where
%%B = 0.
% 
% In the space below, calculate the compressibility factor for a
% van der Waals gas and call it Z_vdw.
% 
>: Z_vdw : Pvdw*V/(R*T)
 
                                     2
                             V (R T V  - a V + a b)
                             ----------------------
                                      3      2
                                R T (V  - b V )
 
% In the space below, calculate the derivative of Z with respect
% to V, and call it dZdV.
% 
>: dZdV : diff(Z_vdw,V)
 
             2
        R T V  - a V + a b   V (2 R T V - a)
        ------------------ + ---------------
               3      2            3      2
         R T (V  - b V )     R T (V  - b V )
                                                 2                2
                                           V (3 V  - 2 b V) (R T V  - a V + a b)
                                         - -------------------------------------
                                                           3      2 2
                                                     R T (V  - b V )
 
% In the space below, calculate the derivative dZ/d(1/V) and call
% it dZd1V. HINT: dy/dx = -(dy/dz)/(dx/dz). 
% 
>: dZd1V : -dZdV/diff(1/V,V)
 
                     2
           2    R T V  - a V + a b   V (2 R T V - a)
        - V  (- ------------------ - ---------------
                       3      2            3      2
                 R T (V  - b V )     R T (V  - b V )
                                                2                2
                                          V (3 V  - 2 b V) (R T V  - a V + a b)
                                        + -------------------------------------)
                                                          3      2 2
                                                    R T (V  - b V )
 
% In the space below, calculate B as the limit of dZ/d(1/V) as V
% approaches infinity. Name it B.
% 
>: B : limit(dZd1V,V,inf)
 
                                     b R T - a
                                   - ---------
                                        R T
 
% In the space below, calculate the Boyle temperature as the temperature
% where B = 0. Call the Boyle temperature TB.
% 
>: TB: rhs(solve(B=0,T)[1])
 
                                        a
                                       ---
                                       b R
 
% Now let's see how well this equation works by looking at trends
% in Boyle temperature for the noble gases. Because we'll be calculating
% numbers and making graphs, we'll switch from Maxima back to Euler.
% 
% In the space below, define arrays a and b that contain the van
% der Waals constants for the noble gases He, Ne, Ar, Kr, and Xe,
% in that order. Take the numbers from Atkins Table 1.6 in the Data
% Section at the end of the book.
% 
>a = [0.0341,0.205,1.337,5.125,4.137]; 
>b = [2.38e-2,1.67e-2,3.20e-2,1.06e-2,5.16e-2];
% In the space below, assign R in units consistent with the units
% for a and b.
% 
>R = 0.082059;
% Next we will calculate an array of predicted Boyle temperatures
% for the noble gases, using the formula we derived in Maxima. Although
% the formula is so simple that it's easy to just retype it in Euler,
% let's see how you would move equations and expressions from Maxima
% into Euler.
% 
% The Euler function mxm() returns a named object in Maxima as a
% text string:
% 
>mxm("TB")
 a/(b*R)
% 
% To evaluate a text string in Euler, use the evaluate() function:
% 
>TB = evaluate(mxm("TB"))
      17.46027991133      149.5929648399      509.1610914098      5891.987058552      977.0338245001 
% 
% In the space below, plot the Boyle temperatures on the y axis
% with atomic number on the x axis. Insert the plot into this notebook.
% As always, title the plot and label the axes.
>atomicnumber = [2,10,18,36,54];
>plot2d(atomicnumber,TB);
>title("Boyle Temperature Trend for the Noble Gases");
>xlabel("Atomic Number");
>ylabel("Boyle Temperature (K)");
>insimg;

% 
% We'd expect the Boyle temperature to increase with increasing
% intermolecular attraction (a) and with decreasing molecular volume
% (b), so it's conceivable that the Boyle temperature might gradually
% rise and then fall going from He to Xe- but the huge spike in
% Boyle temperature at Kr is hard to miss. Either something is wrong
% with the data or calculations, or we have discovered something
% interesting about Kr.
% 
% In the space below, draw the plot again, but this time use van
% der Waals parameters from some other source (such as the CRC handbook).
>a = [0.03412,0.2107,1.345,2.318,4.194];
>b = [2.37e-2,1.709e-2,3.219e-2,3.78e-2,5.105e-2];
>TB = a/(b*R);
>plot2d(atomicnumber,TB);
>title("Boyle Temperature Trend for the Noble Gases");
>xlabel("Atomic Number");
>ylabel("Boyle Temperature (K)");
>insimg;

% 
% The moral of this last plot is: always double-check data points
% that don't fit the pattern implied by the rest of the data. Discordant
% points indicate either an error or an interesting discovery. 
% 
% The little cusp at Ar is interesting. To see if it's real, add
% a curve that shows the true Boyle temperatures as a thicker line
% on the plot. Insert the plot into this notebook again in the space
% below.
% 
>TB = [22.64,122.1,411.5,575,768];
>plot2d(atomicnumber,TB,add=1,thickness=2); insimg;

% 
% Explain the deviation between the true and van der Waals Boyle
% temperatures in terms of intermolecular forces. Does the van der
% Waals equation overestimate or underestimate the strength of intermolecular
% attractions?
% 
% ANSWER: The equation for the Boyle temperature we derived above
% shows that the larger a is, the higher the Boyle temperature will
% be. This suggests that the that the stronger intermolecular attractions
% are, the higher the Boyle temperature must be.
% 
% However the van der Waals predictions of the Boyle temperature
% are too high. That suggests that the van der Waals equation overestimates
% the strength of intermolecular attractions.
% 
—————————-snip here———————————————

Notebook 1.3: Properties of Real Gases from the Ideal Gas Law

Real molecules attract and repel each other. Ideal gas molecules don’t. That’s what causes nearly all deviations from ideality.

The ideal gas law is derived by ignoring attractions and repulsions between molecules. It doesn’t work at all for solids and liquids, which have very strong intermolecular forces. It generally isn’t accurate for real gases, either.

We can still apply the ideal gas law in the limit of low pressure for real gases. At low pressure, the gap between molecules widens, weakening the forces acting between molecules until they are essentially negligible, so the equation applies.

The Euler notebook posted below asks students to connect data for real gases collected at several different pressures with the ideal gas law by extrapolating to zero pressure. It introduces linear regression with Euler.

To use the notebook:

* Download and install the Euler Math Toolbox.
* Cut and paste the code below into a plain text file.
* Save the file with an .en file extension.
* Double-click the file. It should open up in the Toolbox.

I used this in my Fall ‘09 PChem class (it’s the fourth notebook in the series). Italics mark the items students had to fill in themselves. They seemed to have little trouble with this one.

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition.

—————————-snip here———————————————
% Notebook 1.3: Extrapolation to zero pressure
%
% In this notebook, we’ll calculate properties of real gases from the
% ideal gas law by extrapolating to zero pressure.
%
% Consider the following experimental data collected for gaseous
% dimethyl ether at 25 degrees Celsius. The pressures P are in pascals,
% and the densities d are in kg m^-3. R is in Pa m^3/mol K and T is
% K.
%
>P = [12.223,25.20,36.97,60.37,85.23,101.3]*1000;
>d = [0.225,0.456,0.664,1.062,1.468,1.734];
>R = 8.314;
>T = 25 + 273.15;
%
% Let’s find the molar mass of dimethyl ether from this data.
%
% The molar mass M of an ideal gas is mass per mole. The mass of the gas
% is its density times its volume, so
%
% M = dV/n
%
% For an ideal gas, PV=nRT or V/n = RT/P, so
%
% M = dRT/P
%
% The ideal gas law holds when P approaches zero, so for a real gas,
%
% M = lim dRT/P
% P->0
%
% This is called a limiting law. How can we apply the limiting law with the
% pressure and density measurements in P and d?
%
% First, calculate an array of dRT/P values.
%
>dRToverP = d*R*T/P;
%
% Let’s plot these values against P to see if we can linearly extrapolate
% them to zero.
% * When plotting experimental data, it is good practice to show
% the points on the graph (use the points=1 option).
%
>plot2d(P,dRToverP,points=1,style=”[]”);
>title(“Determining the Molar Mass of Dimethyl Ether”);
>xlabel(“Pressure, in Pa”);
>ylabel(“dRT/P, in kg/mol”);
>insimg;

%
% The data looks reasonably linear. It’s worth trying to fit a
% line to it. We can use Euler’s polyfit() function to do this. The
% call looks like
% a = polyfit(x,y,n)
% where:
%
% a is an array of coefficients in the fitted polynomial
% a[1] + a[2]*x + a[3]*x^2 + … + a[n+1]*x^n
% x and y are the arrays of x and y data
% n is the degree of the polynomial to fit (1 for a line)
%
>fit = polyfit(P,dRToverP,1)
0.04587503092209 -3.583072735655e-008
%
% The first number is the intercept; the second is the slope.
%
% To calculate an array of points along the regression line, use the
% polyval() function. The call looks like
% y = polyval(a,x)
% where x, y, and a have been defined above.
%
>x = [0,max(P)];
>y = polyval(fit,x);
%
% Now insert the plot, with the raw data points marked as boxes,
% and the regression line displayed as a solid line.
%
>plot2d(x,y);
>plot2d(P,dRToverP,points=1,style=”[]”,add=1);
>title(“Determining the Molar Mass of Dimethyl Ether”);
>xlabel(“Pressure, in Pa”);
>ylabel(“dRT/P, in kg/mol”);
%
% The y-intercept is M, the value of dRT/P at P=0.
%
>M = fit[1]
0.04587503092209
%
% We can label the y-intercept using the label function. The call
% looks like label(text,x,y,offset). You can play around with the
% offset parameter until the label is where you want it.
%
>label(“M”,0,M,0);
>insimg;

% —————————————————————–
%
% Apply these techniques to solve problem 1.3 in the text, which
% finds absolute zero on the Celsius temperature
% scale from experimental data.
%
% Enter the P and alpha data in the space below.
% (HINT: look carefully at the units.)
%
>P = [749.7, 599.6, 333.1, 98.6];
>alpha = [3.6717, 3.6697, 3.6665,3.6643]*1e-3;
%
% It will be easiest to find alpha for an ideal gas first.
% Do this with the polyfit() function in the space below.
% Call the alpha value for an ideal gas ‘alphaideal’.
>fit = polyfit(P,alpha,1)
0.003662977448051 1.139259281072e-008
>alphaideal = fit[1];
%
% In the space below, insert a plot that extrapolates alpha to zero
% pressure. Mark the data on the line as circles (use style=”o”) and
% add a regression line to the plot. The y-intercept should be shown
% and labeled as “alpha ideal”.
%
>x = [0,max(P)];
>y = polyval(fit,x);
>plot2d(x,y);
>plot2d(P,alpha,points=1,style=”o”,add=1);
>title(“Determining alpha for an ideal gas”);
>xlabel(“Pressure, in torr”);
>ylabel(“alpha, in deg. C^-1″);
>label(” alpha ideal”,20,fit[1],-25);
>insimg;

%
% Now solve Charles’ Law in the form V=Vo(1+alpha*theta) for
% absolute zero. Call the solution ‘absolutezero’.
% HINT: You won’t need an actual value for Vo, if you do your
% algebra correctly.
% HINT: At V=0, theta is absolute zero on the Celsius scale.
% HINT: Charles’ Law applies exactly only to ideal gases.
%
>absolutezero = -1/alphaideal
-273.0019537882
%
% ———————————————————————–
%
% Apply these techniques again to solve problem 1.7(b) in the text,
% which estimates the value of the ideal gas law constant R from
% experimental data.
%
% Enter the P, V, and density(d) data in the space below.
%
>P = [0.750000, 0.500000, 0.250000];
>V = [29.9649, 44.8090, 89.6384];
>d = [1.07144, 0.714110, 0.356975];
%
% Calculate an array of R values calculated at each of the three points.
% Call it RP (for “R at pressure P”).
%
>M = 2*15.9994;
>T = 273.15;
>RP = P*M/(d*T);
%
% In the space below, insert a plot of the RP values as a function
% of pressure. Show the data points as boxes (use style=”[]”). Draw
% the regression line on the plot. AS ALWAYS title the plot and
% label the axes. The axis labels should include the units.
%
>fit = polyfit(P,RP,1)
0.08206187029697 -7.886858741291e-005
>x = [0,max(P)];
>y = polyval(fit,x);
>plot2d(x,y);
>plot2d(P,RP,points=1,style=”[]”,add=1);
>title(“Determination of R”);
>xlabel(“Pressure, in atm”);
>ylabel(“P*M/(density*T), in dm^3 atm/mol K”);
>label(“R”,0,fit[1],0);
>insimg;

%
% In the space below, print your calculated R value.
>R = fit[1]
0.08206187029697

Notebook 1.2 – Solving Equations of State (Molar Volumes of Model Gases)

The Euler notebook posted below asks students to solve various gaseous equations of state for molar volume. This is a common task that’s quite easily done with a computer algebra system, but difficult and tedious to do with pencil and paper.

To use the notebook:

• Download and install the Euler Math Toolbox.
• Cut and paste the code below into a plain text file.
• Save the file with an .en file extension.
• Double-click the file. It should open up in the Toolbox.

I used this in my Fall ‘09 PChem class (the third in the series). Italics mark the items students had to fill in themselves. They had less trouble with this one, although command line syntax was still a source of confusion and frustration. They wanted to type something like root(R*T/V-P,V) or root(RT/V-P,V) instead of root(“R*T/V-P”,V).

All of the notebooks in this series are specifically keyed to Atkins’ Physical Chemistry, 8th edition.

----------------------------snip here---------------------------------------------
Notebook 1.2 - Molar Volume of Model Gases (Solving Equations of State)

%
% In this task, we'll use Euler to numerically solve equations of state.
% Study Example 1.4 on page 18 of Atkins, which estimates the molar
% volume of CO2 at 500 K and 100 atm by treating CO2 as a van der Waals
% gas.
%
% To solve this problem in Euler, set up constants in any convenient
% units (in this case, dm^3 for volume, atm for pressure, and K for
% temperature:
>R = 0.082059;
>a = 3.592;
>b = 4.267e-2;
>T = 500;
>P = 100;
%
% Now cast the equation into a form that looks like f(x) = 0.
% Atkins rearranges the van der Waals equation into polynomial
% in V:
% V^3-(b+R*T/P)*V^2 + (a/P)*V - a*b/P = 0.
% but that really isn't necessary. *Any* rearrangement of the equation
% we want to solve that puts zero on one side of the equation will work.
%
% To solve the equation in Euler, we'll need an initial guess for V.
% The ideal gas volume will do.
>V = R*T/P
0.410295
%
% Now use Euler's 'root' function, which looks like
% root("expression",x)
% where expression (enclosed in double quotes) is the nonzero side of the
% equation we want to solve, and x is the variable in that expression we
% want to solve for.
%
>root("V^3-(b+R*T/P)*V^2 + (a/P)*V - a*b/P",V)
0.3663332616871
%
% Note that root has reset the value of V:
%
>V
0.3663332616871
%
% Atkins started with the equation
% P = R*T/(V-b) - a/V^2
% We could have subtracted P from both sides of this equation to more
% simply cast it into the f(x) = 0 form:
%
>root("R*T/(V-b)-a/V^2-P",V)
0.3663332616871
%
% Atkins claims that a perfect gas under these conditions has a molar
% volume of 0.410 dm^3/mol. We can verify this by typing
%
>root("R*T/V - P", V)
0.4102949999999
% ---------------------------------------------------------------------
% TASK 3. Use this second approach to calculate the molar volume of
% argon in dm^3 per mole at 100 degrees Celsius and 100 atm, on the
% assumption that it's
% A) an ideal gas;
% B) a hard sphere gas (assuming that each molecule is a
% sphere with radius of 0.15 nm);
% C) a van der Waals gas;
% D) a Dieterici gas;
% E) a virial gas (ignoring all virial coefficients beyond B).
%
>T = 100+273.15;
>P = 100;
%
% 3A) The molar volume as an ideal gas:
%
>V = R*T/P
0.3062031585
%
% 3B) The molar volume as a hard sphere gas, with radius 0.15 nm:
% HINT: In Euler, the constant pi is written as %pi
% HINT: What units should you put b into
% if V is in dm^3/mol?
%
>r = 0.15*1e-9/1e-1;
>NA = 6.022e23;
>b = 4*(4/3)*%pi*r^3*NA;
>root("R*T/(V-b)-P",V)
0.3402567662277
%
% 3C) The molar volume as a van der Waals gas:
%
>a = 1.337;
>b = 3.20e-2;
>root("(P+a/V^2)*(V-b)-R*T",V)
0.2981759820445
%
% 3D) The molar volume as a Dieterici gas:
% HINT: See Table 1.7 in Atkins.
% HINT: a and b for the Dieterici gas are NOT the same as
% a and b for a van der Waals gas!
% HINT: All of the expressions in the Critical Constants columns
% in Table 1.7 are equal.
% HINT: In Euler, the constant e is written as %e
% HINT: In Euler, e^x can be written as exp(x)
%
>B = 3*b/2;
>A = 4*%e^2*B^2*a/(27*b^2);
>root("R*T*exp(-A/(R*T*V))/(V-B)-P",V)
0.245635997406
%
% 3E) The molar volume as a virial gas (with all virial coefficients
% beyond B ignored)
% HINT: Watch your units; what should the units of B be, if V
% should be in dm^3/mol?
% HINT: B varies with temperature.
%
>B = -4.2/1000;
>root("(R*T/V)*(1+B*P)-P",V)
0.17759783193

----------------------------snip here---------------------------------------------

Filling in HTML attributes with XSL

Problem: I want to transform some XML that looks like

<flash href="drill-1-scientific-notation" width="550" height="400"/>

into browser-independent HTML, using XSL. The output code has to copy the XML attributes into HTML attributes. I can’t use xsl:value-of tags to do that, because I can’t put tags inside of tags.

Solution: use {@attribute} wherever the value of an attribute has to be inserted. The XSL template for my flash tag looks like

<xsl:template match="flash">
<object classid="clsid:d27cdb6e-ae6d-11cf-96b8-444553540000" codebase="http://download.macromedia.com/pub/shockwave/cabs/flash/swflash.cab#version=9,0,0,0" width="{@width}" height="{@height}" id="{@href}" align="middle">
<param name="allowScriptAccess" value="sameDomain" />
<param name="allowFullScreen" value="false" />
<param name="movie" value="{@href}.swf" />
<param name="quality" value="high" />
<param name="bgcolor" value="#ffffff" />
<embed src="{@href}.swf" quality="high" bgcolor="#ffffff" width="{@width}" height="{@height}" name="{@href}" align="middle" allowScriptAccess="sameDomain" allowFullScreen="false" type="application/x-shockwave-flash" pluginspage="http://www.macromedia.com/go/getflashplayer" />
</object>
</xsl:template>