* Problem 6.30. Para and ortho hydrogen calculation. * Pat Arnott. November 8, 2011. IMPLICIT DOUBLE PRECISION (a-h,o-z) ! my crutch for defining variables in the traditional way. DOUBLE PRECISION k,kT,T(200000),lnZeven(200000),lnZodd(200000),beta(200000),Ceven(200000),Codd(200000) DOUBLE PRECISION Zeven(200000),Zodd(200000),Eeven(200000),Eodd(200000) eps=0.0076d0 ! eV k=8.617d-5 ! eV/K Boltzmann's constant. imax=10000 ! This is a guess at the resolution needed to give a good finite difference derivative. DO i=1,imax T(i) = 2.0d0 + 490.0d0*dFLOAT(i)/dFLOAT(imax) ! Temperature in Kelvin. kT=k*T(i) beta(i)=1.0d0/kT Zeven(i)=0.0d0 ; Zodd(i)=0.0d0 DO j=0,100,2 ! Even terms sum. xj=dFLOAT(j) Zeven(i) = Zeven(i) + (2.0d0*xj+1.0d0)*dEXP(-xj*(xj+1.0d0)*eps/kT) ENDDO DO j=1,101,2 ! Odd terms sum. xj=FLOAT(j) Zodd(i) = Zodd(i) + (2.0d0*xj+1.0d0)*dEXP(-xj*(xj+1.0d0)*eps/kT) ENDDO lnZeven(i)=dLOG(Zeven(i)) lnZodd(i)= dLOG(Zodd(i)) ENDDO * Calculate E and then C. E=-PartialLNz / PartialBeta. C=PartialE / PartialT DO i=2,imax-1 Eeven(i)=-Derivative(lnZeven(i+1),lnZeven(i),lnZeven(i-1),beta(i+1),beta(i),beta(i-1)) Eodd(i)=-Derivative(lnZodd(i+1),lnZodd(i),lnZodd(i-1),beta(i+1),beta(i),beta(i-1)) ENDDO Eeven(1)=Eeven(2) ; Eodd(1)=Eodd(2) ; Eeven(imax)=Eeven(imax-1) ; Eodd(imax)=Eodd(imax-1) ; DO i=2,imax-1 Ceven(i)=Derivative(Eeven(i+1),Eeven(i),Eeven(i-1),T(i+1),T(i),T(i-1)) Codd(i)=Derivative(Eodd(i+1),Eodd(i),Eodd(i-1),T(i+1),T(i),T(i-1)) ENDDO Ceven(1)=Ceven(2) ; Codd(1)=Codd(2) ; Ceven(imax)=Ceven(imax-1) ; Codd(imax)=Codd(imax-1) ; * Calculate average values for C from these high resolution calculations. * Then print out the results. OPEN(1,FILE='problem6_30.txt') WRITE(1,*) 'kT/eps Cpara/k Cortho/k Chydrogen/k' DO i=1,imax,2*imax/100 Tave=0.0d0 ; Ce=0.0d0 ; Co=0.0d0 ; iii=0d0 DO ii=0,2*imax/100-1 iii=iii+1 Tave=Tave+T(i+ii) ; Ce=Ce+Ceven(i+ii) ; Co=Co+Codd(i+ii) ENDDO Tave=Tave/dFLOAT(iii) ; Ce=Ce/dFLOAT(iii) ; Co=Co/dFLOAT(iii) x=k*Tave/eps WRITE(1,*) x,Ce/k,Co/k,(0.25d0*Ce+0.75d0*Co)/k ENDDO CLOSE(1) END Function Derivative(yip1,yi,yim1,xip1,xi,xim1) IMPLICIT Double Precision (a-h,o-z) Term1=(yip1-yi)/(xip1-xi) ; Term2=(yi-yim1)/(xi-xim1) Derivative=(Term1+Term2)/2.0d0 RETURN END