Assignment:
Write a
code to evaluate for the following simple model of matter.
Two and 3d integration
The integral of interest is
This defines
Where f is a simple approximation to the density distribution in a simple cubic crystal given by
The lowest level function is
FUNCTION FUNP(PHI) for\FUNP.FOR
IMPLICIT REAL*8 (A-H,O-Z)
COMMON /PASS/R,THETA
Z=R*COS(THETA)
RST=R*SIN(THETA)
X=RST*COS(PHI)
Y=RST*SIN(PHI)
FUNP=FUN(X,Y,Z)
RETURN
END
FUNCTION FUN(X,Y,Z)
DATA PI/3.141592653589790D0/
DATA A,B,C/1D0,1D0,1D0/
IF(ABS(X).LT..5D0*A.AND.ABS(Y).LT..5D0*B.
2 AND.ABS(Z).LT..5D0*C)THEN
FUN=0
ELSE
2 *COS(PI*C*Z)
FUN=8*
ENDIF
RETURN
END
The C version of this is cpp\funp.cpp
double fun(double x,double y,double z)
{double a=1,b=1,c=1,pi=3.141592653589790,
retval=0;
if((abs(x) < .5*a) && (abs(y) < .5*b) && (abs(z) < .5*c))return retval;
retval=cos(pi*x)*cos(pi*y)*cos(pi*z);
retval*=retval*8;
return retval;}
double funp(double phi)
{double x,y,z,rst,retval;
z=pr*cos(ptheta);
rst=pr*sin(ptheta);
x=rst*cos(phi);
y=rst*sin(phi);
retval=fun(x,y,z);
return retval;}
A good deal
of physics can be added to this by using random numbers to vary the locations
of the individual locations of the atoms.
This can even be iterated on with various temperatures. In all honesty, however, the notion that r
dominates is usually used to find by averaging over the r's found at different
distances.
To be completely specific, the next higher level function is for\FUNT.FOR
FUNCTION FUNT(THETA)
IMPLICIT REAL*8 (A-H,O-Z)
COMMON /PASS/R,PTHETA
DATA PI/3.141592653589790D0/
PTHETA=THETA
M=100
FUNT=0
H=2*PI/M
DO I=1,M
PHI=H*(I-.5D0)
FUNT=FUNT+FUNP(PHI)
ENDDO
FUNT=SIN(THETA)*H*FUNT
RETURN
END
Note that the theta passed to the lower level function needs to be given a different name in order to go from the argument list to the common list. The C version of this is cpp\funt.cpp
double funt(double theta)
{int m=100,i;
double h,pi=3.141592653589790,
retval,phi;
ptheta=theta;
h=2*pi/m;
retval=0;
for (i=0;i<m;++i)
{phi=h*(i+.5);
retval+=funp(phi);}
retval*=sin(theta)*h;
return retval;}
In C the sketch for a code is in cpp\main.cpp
#include <stdio>
#include <math.h>
double pr=5.,ptheta;
#include "funp.cpp"
#include "funt.cpp"
void main(void)
{double theta=.5,ft;
FILE *rout;
if((rout=fopen("test.out","w"))==NULL){printf("open err test.out\n");}
ft=funt(theta);
fprintf(rout," %20.16lf %20.16lf\n",theta,ft);
return;}
Note that the external variable pr and ptheta substitute for the common /pass/ in the Fortran. These codes are meant as a starting point only. In this assignment, you will need to write code. You will need to test the functions. For example for\TFUNP.FOR
IMPLICIT REAL*8 (A-H,O-Z)
COMMON /PASS/R,THETA
DATA PI/3.141592653589790D0/
R=3.35
THETA=1.32
H=2*PI/127
OPEN(1,FILE='FUNP.OUT')
DO I=1,127
PHI=(I-.5D0)*H
FT=FUNP(PHI)
WRITE(1,*)PHI,FT
ENDDO
END
C$INCLUDE FUNP
Figure 1 Lowest level integrand for r=3.35, theta = 1.32, versus phi.
Documented code for calculating ρ(r) in C with a make file is in cpp\twobody.tar. The output, included as output.txt is plotted below.
Figure 2 The function rho(r) versus r for this model
The Laurent transform
..\..\wsteve\fmonte2\FMONTE.HTM
The Integral of the Lagrange polynomial through the points a, b, c and d.
Consider the Fourier transform from begr to endr of a function defined as a Lagrange polynomial.
The function through these points is
The integral of f is
The Integrals are of the form
Define and let
so that the integral becomes
Define
So that
The integral is
The code to integrate from the first point to the last in Fortran is in for\AINT.FR
FUNCTION TAINT(XDAT,FDAT,N)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION XDAT(*)
COMPLEX*16 FDAT(*)
C EVALUATES THE INTEGRAL FROM XDAT(1) TO
XDAT(N)
C LAGRANGE POLYNOMIAL THROUGH THE
NEAREST 4 POINTS
BEGR=XDAT(1)
ENDR=XDAT(3)
IF(N.LE.4)THEN
IF(N.LT.4)THEN
PRINT*,' AT LEAST 4 POINTS ARE NEEDED
FOR'
PRINT*,' INTERPOLATION'
READ(*,*)ITEST
STOP
ELSE
ENDR=XDAT(4)
ENDIF
ENDIF
TAINT=AINT2(XDAT(1),XDAT(2),XDAT(3),
A XDAT(4),BEGR,ENDR
2 ,FDAT(1),FDAT(2),FDAT(3),FDAT(4))
IF(N.EQ.4)RETURN
DO I=3,N-3
BEGR=XDAT(I)
ENDR=XDAT(I+1)
TAINT=TAINT+
2 AINT2(XDAT(I-1),XDAT(I),XDAT(I+1),XDAT(I+2),BEGR,ENDR,
3 FDAT(I-1),FDAT(I),FDAT(I+1),FDAT(I+2))
ENDDO
BEGR=XDAT(N-2)
ENDR=XDAT(N)
TAINT=TAINT+
2 AINT2(XDAT(N-3),XDAT(N-2),XDAT(N-1),XDAT(N),BEGR,ENDR,
3 FDAT(N-3),FDAT(N-2),FDAT(N-1),FDAT(N))
RETURN
END
FUNCTION AINT2(SA,SB,SC,SD,BEGR,ENDR,FA,FB,FC,FD)
IMPLICIT REAL*8 (A-H,O-Z)
COMPLEX*16 FA,FB,FC,FD
H=ENDR-BEGR
XMID=(BEGR+ENDR)/2
AINT2=FA*AIHAT2(SB,SC,SD,H,XMID)/((SA-SB)*(SA-SC)*(SA-SD))
2 +FB*AIHAT2(SA,SC,SD,H,XMID)/((SB-SA)*(SB-SC)*(SB-SD))
3 +FC*AIHAT2(SA,SB,SD,H,XMID)/((SC-SA)*(SC-SB)*(SC-SD))
4 +FD*AIHAT2(SA,SB,SC,H,XMID)/((SD-SA)*(SD-SB)*(SD-SC))
RETURN
END
C
FUNCTION AIHAT2(SB,SC,SD,H,XMID)
IMPLICIT REAL*8 (A-H,O-Z)
CAPB=(XMID-SB)
CAPC=(XMID-SC)
CAPD=(XMID-SD)
AIHAT2=H*(CAPB*CAPC*CAPD+H*H*
2(CAPB+CAPC+CAPD)/12)
RETURN
END
Any function can be found using findfun and integrated by the above to yield
In the above let so that
The function G is known on the grid of points found by findfun. In a general integral of the form
make the variable change using a positive definite g
Then