-
Notifications
You must be signed in to change notification settings - Fork 0
/
ideal-eutectic
51 lines (51 loc) · 2.37 KB
/
ideal-eutectic
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#!/usr/bin/perl
&setup;$file_count=@ARGV;$terminate=1 if $file_count<2;$file_count--;
foreach $file(@ARGV){
if (-e $file){
push @mols,$file;
$head=`head -1 $file`;
($temperature,$enthalpyfus)=grep length, split(/\s+/,$head);
$terminate++ if $temperature<0; #melt temp must be in K, so positive
$terminate++ if $enthalpyfus<0; #enthalpy of fusion is positive for everything but helium...
#print "$file $temp $enthfus\n";
push @temps,$temperature;push @enthalpies,$enthalpyfus;
}else{
$terminate++; warn "cannot find file: $file";
};
};
die "expecting input of form:\n./exe [molecule 1] [molecule 2] [opt more molecules (as many as desired)]\nnote that T_melt and dH_fusion are expected to be positive" unless ($terminate==0);
foreach $index (0..$file_count-1){ @concs[$index]=1/($file_count+1); };
while($diff>$tol){
print "iteration: $iter\n" if $debug==1;
$index=0;$eutectic_est=1/(1/@temps[$index]-$gas_const*log(@concs[$index])/@enthalpies[$index]);
print "$iter eutectic est $index @concs[$index] $eutectic_est\n" if $debug==1;
$aoo=0;
foreach $index(1..$file_count){
@concs[$index]=exp(-@enthalpies[$index]/$gas_const*(1/$eutectic_est-1/@temps[$index]));
$aoo+=@concs[$index]; #add up the other mole fractions
print "$index mol frac: @concs[$index] sum: $aoo\n" if $debug==1;
};
$aoo=1-$aoo; #subtract from 1 to get new estimate for molecule 1 mole fract
$foo=(@concs[0]+$aoo)/2; #average
print "the new est. mol. fract. for @mols[0] is $foo from $aoo\n" if $debug==1;
$diff=abs($foo-$aoo)/($foo)*100;
print "$iter gen diff $diff\n" if $debug==1;
$index=0;@concs[$index]=$foo;
$iter++;
die "maximum iterations reached: $max_iter\n" if $iter>$max_iter;
};
print "********\neutectic summary: $eutectic_est K\n";
foreach $index (0..$file_count){
print "@mols[$index] mole fraction: @concs[$index]\n";
};
print "convergence tolerance in mole fraction: $tol
required $iter iterations to converge within $diff\n********\n";
sub setup {
$debug=0;
$terminate=0;
$gas_const=0.0083144598; #kJ/mol/K
$tol=0.00001;
$diff=100;
$iter=0;
$max_iter=500;
};