#!/usr/bin/perl

#Problem 6.30.  A Perl script to calculate para and ortho hydrogen.
#Written by Ben Sumlin, November 2011.
#PHYS425 Homework 14

open (FILE, '>>problem630.txt');
print FILE "kT/eps Cpara/k Cortho/k Chydrogen/k \n";

#Declare some constants...
$eps=0.0076; # eV
$k=0.00008617; # Boltzmann's constant in eV/K

#Declare some arrays.  Perl doesn't care about array size at initialization.
@T;
@Beta;
@lnZeven;
@lnZodd;
@beta;
@Ceven;
@Codd;
@Zeven;
@Zodd;
@Eeven;
@Eodd;

$imax=10000; #resolution adjustment

#The foreach command is more elegant to deal with arrays in Perl, but since
#we didn't define the size of arrays at initialization, we have to use the for command.

for($i=0; $i<$imax; $i++)
{
	$T[$i]=2.0 + 490.0*($i/$imax);
	$KT=$k*$T[$i];
	$Beta[$i]=1/$KT;
	
	$Zeven[$i]=0;
	$Zodd[$i]=0;

	for($j=0; $j<100; $j=$j+2) #The even values
	{
		$Zeven[$i]=$Zeven[$i]+(2*$j+1) * exp(-$j*($j+1)*$eps/$KT);
	}

	for($j=1; $j<101; $j=$j+2) #The odd values
	{
		$Zodd[$i]=$Zodd[$i]+(2*$j+1) * exp(-$j*($j+1)*$eps/$KT);
	}
	
	$lnZeven[$i]=log($Zeven[$i]);
	$lnZodd[$i]=log($Zodd[$i]);
}

for($i=1; $i<$imax-1; $i++) #Calculate <E>
{
	$Eeven[$i]=-((($lnZeven[$i+1]-$lnZeven[$i])/($Beta[$i+1]-$Beta[$i]))+(($lnZeven[$i]-$lnZeven[$i-1])/($Beta[$i]-$Beta[$i-1])))/2;
	$Eodd[$i]=-((($lnZodd[$i+1]-$lnZodd[$i])/($Beta[$i+1]-$Beta[$i]))+(($lnZodd[$i]-$lnZodd[$i-1])/($Beta[$i]-$Beta[$i-1])))/2;
}
$Eeven[0]=$Eeven[1];
$Eodd[0]=$Eodd[1];
$Eeven[$imax]=$Eeven[$imax-1];
$Eodd[$imax]=$Eodd[$imax-1];

for($i=1; $i<$imax-1; $i++) #Calculate C
{
	$Ceven[$i]=((($Eeven[$i+1]-$Eeven[$i])/($T[$i+1]-$T[$i]))+(($Eeven[$i]-$Eeven[$i-1])/($T[$i]-$T[$i-1])))/2;
	$Codd[$i]=((($Eodd[$i+1]-$Eodd[$i])/($T[$i+1]-$T[$i]))+(($Eodd[$i]-$Eodd[$i-1])/($T[$i]-$T[$i-1])))/2;
}

$Ceven[0]=$Ceven[1];
$Codd[0]=$Codd[1];
$Ceven[$imax]=$Ceven[$imax-1];
$Codd[$imax]=$Codd[$imax-2];

for($i=0; $i<$imax; $i=$i+2*$imax/100)
{
	$Tave=0;
	$Ce=0;
	$Co=0;
	$iii=0;
	for($ii=0; $ii<2*$imax/100 - 1; $ii++)
	{
		$iii=$iii+1;
		$Tave=$Tave+$T[$i+$ii];
		$Ce=$Ce+$Ceven[$i+$ii];
		$Co=$Co+$Codd[$i+$ii];
	}
	$Tave=$Tave/$iii;
	$Ce=$Ce/$iii;
	$Co=$Co/$iii;
	$x=$k * $Tave/$eps;
	$y=$Ce/$k;
	$z=$Co/$k;
	$rawr=(0.25*$Ce+0.75*$Co)/$k;

	print FILE $x;
	print FILE " ";
	print FILE $y;
	print FILE " ";
	print FILE $z;
	print FILE " ";
	print FILE $rawr;
	print FILE "\n";
}

close (FILE);