kill(all)$

numer:true;

maxima:lab:vapor pressure ,version 1.0
peter.vlasschaert@gmail.com,27/11/2019

load : packages
=============

load(lrats)$

definition : vapor pressure
====================
log10(P0) = A+B/(t+230)
------------------------------------
A,B = constant
  t    = temperature

qq:A+B/(t+230);

data: genmatrix(a, 9, 6);

data[1,1]:xx;

data;

How to insert : Data , computation into matrix data.
========================================

three columns : data    =>  list : 'form'
------------------------------------------------

definition : log with base : 10
----------------------------------------

log10(x):=float(log(x)/log(10));

check definition
============

log10(100);

definition : vapor pressure
------------------------------------
B : negative, 'see above'

 vap(A,B,t):=''qq;

data : methanol : CH3OH
-----------------------------------

info :
===

dd:["no","t(°C)","P0(mmHg)","log10(P0)","b*10^3","l"];

data:
====

nd:[5,6,7,9,10,12,13,15];
Cd:[49.4,45.9,43.1,32.1,26.3,18.65,15.0,8.0];
Pd:[399.3,345.5,303.3,180.4,134.9,90.6,73.7,49.5];

for i:1 thru length(dd) do(
data[1,i]:dd[i]
   )$

for i:2 thru length(nd)+1 do(
data[i,1]:nd[i-1],
data[i,2]:Cd[i-1] ,
data[i,3]:Pd[i-1],
data[i,4]:log10(Pd[i-1])       
   )$   

vapor pressure : CH3OH
====================

data;

y: values
=======
y = log10(P0)

p1:map(log10,Pd);

x:values
======
x=1/(t+230)

p2:makelist(subst([t=Cd[i]],qq),i,1,length(Cd));

p2a:coeff(p2,B);

list of data points
==============

xy:makelist([p2a[i],p1[i]],i,1,length(p1));

relation between :
=============
1e)  log10(P0) = A+B/(t+230)
2e)  y = A+B*x

wxplot2d([['discrete, xy]])$

wxplot2d([['discrete, xy],['discrete, xy]],[style,points,lines]);

help : ? lsquares_estimates;

How to use :lsquares_estimates;

1e) load packages
-------------------------

load(lsquares)$

2e)  convert list to matrix
----------------------------------

p3:matrix(p2a,p1)$

p3a:transpose(p3);

3e) lsquares_estimates
---------------------------------
syntax: lsquares_estimates(p3a,var'list',expression,values 'list')

fpprec:8$

p4:lsquares_estimates(p3a,[x,y],y=A0+B0*x,[A0,B0])$

get equation :  y=A0+B0*x

p4a:float(p4),flatten;

p4b:flatten(p4a);

p4a1:p4b[1];
p4a2:p4b[2];

How to find : b,l
============

log10(P0) = A+B/(t+230)

1e)   How to find the correct value A,B
------------------------------------------------------
ΔA,ΔB  ERROR

p5a: A=A0+ΔA;
p5b: B=B0+ΔB;

2e)  subst qq 1e)
------------------------

p6:subst([p5a,p5b],qq)-y;

p6a:expand(p6);

3e)  find coef ΔA,ΔB,l
=================
y = log10(P0)

p8a2:a=ratcoeff(p6,ΔA);
p8b2:b=ratcoeff(p6,ΔB);
p8c2:l=ev(p6,ΔA=0,ΔB=0);

theory of probability
::::::::::::::::::::::::::::::::::

p7:v=p6a;

p7a:[v=v[i],t=t[i]];

p7b:[a=a[i],b=b[i],l=l[i],t=t[i],v=v[i]];

sum ((v[i])^2) = min , i=1,..,n , n= number measurement

4e)  How to find min from above theory
==============================

p8:v=a*ΔA+b*ΔB+l;

p8a :expand(p8*p8);

summation : with use  f1
---------------------------------

f1(q1,q2):=block([],
    if q1=q2 then 
     ( m1:sum(q1[i]^2,i,1,n),
      return(m1))
    else
      ( m2:sum(q1[i]*q2[i],i,1,n),
       return(m2)));

apply :f1 to 'p8a'
=============

p8b:v^2=f1(b,b)*(ΔB)^2+2*f1(a,b)*ΔA*ΔB+2*f1(b,l)*ΔB+f1(a,a)*(ΔA)^2+2*f1(a,l)*ΔA+f1(l,l);

p8c:lhs(p8b) =subst([a[i]=1],rhs(p8b));

p8c1:rhs(p8c);

rhs 'p8c' :function (ΔA,ΔB)
=====================
theory probability :v^2=min (sum(v[i]^2,i,1,n))

1e)  How to find min of function
=========================

d(f)/dΔA = 0  & d(f)/dΔB = 0

f: p8c1;

p9a:'diff(f,ΔA) = diff(f,ΔA);
p9b:'diff(f,ΔB) = diff(f,ΔB);
p9a1:rhs(p9a) = 0;
p9b1:rhs(p9b) = 0;

HOW TO SOLVE : 'p9a1,p9b1'
========================
solve : 'ΔA,ΔB'

p9c:[ΔA,ΔB];

p9c1:solve([p9a1,p9b1],p9c);

declare(i,integer);

p9c2:flatten(p9c1);

p9c3:subst([n=8],p9c2);

find : b[i],l[i]
=========
you need : lratsubst  & not  subst

ratprint:false$

p10a:p9c3[1];
p10b:rhs(p10a);

p10c:changevar(%, j=i, j, i);

p10d:p10c,simpsum;

p11a:p8b2;

p11b:makelist(b[k]=1/(Cd[k]+230),k,1,8);

p12a:p8c2;

p12b:subst([l=l[k],t=t[k],y=y[k]],p12a);

rem : not use index for list already declared before.
Pd,Cd......

p12c:makelist(ev(p12b),k,1,8);

a)  without  function

p12d:makelist(t[k]=Cd[k],k,1,8);

p12e:subst(p12d,p12c);

p12e1:subst(p4b,p12e);

b) with  function , use 'ev' both side

p12e2:makelist(ev(y[k])=ev(log10(Pd[k])),k,1,8);

p12e3:subst(p12e2,p12e1);

find  ΔA
======
use : p11b,p12e3

p13a:ΔA = p10d;
p13b:subst(p11b,p13a);
p13c:subst(p12e3,p13b);

find  ΔB
======
use : p11b,p12e3

p14a:p9c3[2];
p14b:rhs(p14a),simpsum;

p15a:ΔB = p14b;
p15b:subst(p11b,p15a);
p15c:subst(p12e3,p15b);

data table 
---------------
use :p11b,,p12e3

for i:2 thru length(nd)+1 do(
    data[i,5]:rhs(p11b[i-1])*1000,
    data[i,6]:rhs(p12e3[i-1])    
    )$

p16:data;

correct value:
A = A0+ΔA
B = B0+ΔB
----------------
error =ΔA,ΔB

p17a:A=A0+ΔA;

p7b:p4a1;
p17c:p13c;

p18a:B=B0+ΔB;

p7b:p4a2;
p17c:p15c;

ref : vapour liquid equilibrium(second edition,eduard hala)