program jam3
use physics
use utility
implicit none
integer :: i, j, n
integer, parameter ::
yd=50, xd=50, zd=50
real, parameter :: qe=1.602e-19
real :: q=1.0e10*qe
type (vec3d) :: point1,
delta= vec3d(2,2,2)
character*8 :: filename='evec'
character*40 :: prompt
= 'enter xi, yi, and zi'
real, dimension(xd) ::
x
real, dimension(yd) ::
y
real, dimension(zd) ::
z
real, dimension(xd,yd,zd)
:: v, vsum
type (vec3d), dimension(xd,yd,zd)
:: f, fsum
open(unit=7,file=filename)
print *, prompt
read (*,*) point1%x,point1%y,
point1%z
do i=1,xd
x(i)=delta%x*(i-1)
end do
do j=1,yd
y(j)=delta%y*(j-1)
end do
do n=1,zd
z(n)=delta%z*(n-1)
end do
call electricfield(f,
point1, x, y, z, q)
call potential(v, point1,
x, y, z, q)
fsum=f
vsum=v
point1%x=point1%x+20
point1%z=point1%z-20
call electricfield(f,
point1, x, y, z, -q)
call potential(v, point1,
x, y, z, -q)
vsum=vsum+v
call add(fsum,f)
call write_dx(x,y,z,filename,7)
do n=1, zd
do j=1,
yd
do i=1, xd
write (7,*) fsum(i,j,n)
write (8,*) vsum(i,j,n)
end do
end
do
end do
end program jam3
module utility
! this module contains subroutines for writing header files
contains
subroutine write_asf(col,row,unitnum)
! this subroutine writes ASCII special file headers
! col contains column values and row contains row values
! unitnum contains fortran unit number for output
real, dimension(:) :: col,row
integer :: unitnum
write (unit=unitnum,fmt=*) size(col),
size(row)
write (unit=unitnum,fmt=*) 0,0
write (unit=unitnum,fmt=*) row
write (unit=unitnum,fmt=*) col
end subroutine write_asf
subroutine write_dx(x,y,z,name,unitnum)
! header file for IBM data explorer
real, dimension (:) :: x,y,z
integer :: unitnum
character*(*) name
write (unit=unitnum,fmt=*) '# Data Explorer
Header File'
write (unit=unitnum,fmt=*) '#'
write (unit=unitnum,fmt=*) 'object "data"
array items',&
&size(x)*size(y)*size(z)
write (unit=unitnum,fmt=*) ' rank 1
shape 3'
write (unit=unitnum,fmt=*) ' type float'
write (unit=unitnum,fmt=*) ' data file
',trim(name), ',0'
write (unit=unitnum,fmt=*) ' attribute
"dep" string "positions"'
write (unit=unitnum,fmt=*) 'object "pos"
gridpositions counts',&
& size(x), size(y), size(z)
write (unit=unitnum,fmt=*) ' origin
', z(1), y(1), x(1)
write (unit=unitnum,fmt=*) ' delta ',
0, 0, x(2)-x(1)
write (unit=unitnum,fmt=*) ' delta ',
0, y(2)-y(1), 0
write (unit=unitnum,fmt=*) ' delta ',
z(2)-z(1), 0, 0
write (unit=unitnum,fmt=*) 'object "con"
gridconnections counts',&
& size(x), size(y), size(z)
write (unit=unitnum,fmt=*) ' attribute
"element type" string "cubes"'
write (unit=unitnum,fmt=*) ' attribute
"ref" string "positions"'
write (unit=unitnum,fmt=*) 'object "f1"
field'
write (unit=unitnum,fmt=*) ' component
"data" "data"'
write (unit=unitnum,fmt=*) ' component
"positions" "pos"'
write (unit=unitnum,fmt=*) ' component
"connections" "con"'
end subroutine write_dx
end module utility
module physics
! this module contains physics subroutines
type vec2d
real :: x, y
end type vec2d
type vec3d
real :: x, y, z
end type vec3d
! global constants for subroutines
real :: k=8.99E9, rad=1
interface potential
module procedure potential2D,
potential3D
end interface
interface electricfield
module procedure electricfield2D,
electricfield3D
end interface
interface add
module procedure addvec2d,
addvec3d
end interface
contains
subroutine addvec2d(a, b)
! this subroutine adds two 2d vectors
! input: a,b
! output: a
type (vec2d), dimension(:,:) :: a, b
a%x=a%x+b%x
a%y=a%y+b%y
end subroutine addvec2d
subroutine addvec3d(a, b)
! this subroutine adds two 3d vectors
! input: a,b
! output: a
type (vec3d), dimension(:,:,:) :: a,b
a%x=a%x+b%x
a%y=a%y+b%y
a%z=a%z+b%z
end subroutine addvec3d
subroutine electricfield2D(e,point0,x,y,q)
! this subroutine calculates 2D electric field of point charge
! input variables are: point0,q,x,y
! output variables are: e
! point0%x and point0%y are locations of point charge
! q is the charge in coulombs
! x and y are locations in space where electric field is calculated
! e%x and e%y are the components of the electric field
real :: v, r
real :: q
type (vec2d) :: point0
real, dimension(:) :: x, y
type (vec2d), dimension(:,:) :: e
do i=1, size(x)
do j=1, size(y)
r=sqrt((x(i)-point0%x)**2+(y(j)-point0%y)**2)
If (r<rad) r=rad
v=(k*q)/r**2
e(i,j)%x=v*(x(i)-point0%x)/r
e(i,j)%y=v*(y(j)-point0%y)/r
end do
end do
end subroutine electricfield2D
subroutine electricfield3D(e,point0,x,y,z,q)
! this subroutine calculates 3D electric field of point charge
! input variables are: point0,q,x,y,z
! output variables are: e
! Point0%x, point0%y, and point0%z are locations of point charge,
! q is the charge in coulombs
! x, y, and z are locations in space where electric field is calculated
! e%x, e%y, and e%z are the components of the electric field
real :: v, r
real :: q
type (vec3d) :: point0
real, dimension(:) :: x, y, z
type (vec3d), dimension(:,:,:) :: e
do n=1, size(z)
do j=1, size(y)
do i=1, size(x)
r=sqrt((x(i)-point0%x)**2+(y(j)-point0%y)**2+(z(n)-point0%z)**2)
If(r<rad) r=rad
v=(k*q)/r**2
e(i,j,n)%x=v*(x(i)-point0%x)/r
e(i,j,n)%y=v*(y(j)-point0%y)/r
e(i,j,n)%z=v*(z(n)-point0%z)/r
end do
end do
end do
end subroutine electricfield3D
subroutine potential2D(pot,point0,x,y,q)
! this subroutine calculates 2D electrical potential of point charge
! input variables are: point0,q,x,y
! output variables are: pot
! point0%x and point0%y are locations of point charge,
! q is the charge in coulombs
! x and y are locations in space where potential is calculated
! pot is the electrical potential
real :: r
real :: q
type (vec2d) :: point0
real, dimension(:) :: x, y
real, dimension(:,:) :: pot
do i=1, size(x)
do j=1, size(y)
r=sqrt((x(i)-point0%x)**2+(y(j)-point0%y)**2)
If (r<rad) r=rad
pot(i,j)=(k*q)/r
end do
end do
end subroutine potential2D
subroutine potential3D(pot,point0,x,y,z,q)
! this subroutine calculates 3D electrical potential of point charge
! input variables are: point0,q,x,y,z
! output variables are: pot
! point0%x, point0%y, and point0%z are locations of point charge,
! q is the charge in coulombs
! x, y, and z are locations in space where potential is calculated
! pot is the electrical potential
real :: r
real :: q
type (vec3d) :: point0
real, dimension(:) :: x, y, z
real, dimension(:,:,:) :: pot
do n=1, size(z)
do j=1, size(y)
do i=1, size(x)
r=sqrt((x(i)-point0%x)**2+(y(j)-point0%y)**2+(z(n)-point0%z)**2)
If (r<rad) r=rad
pot(i,j,n)=(k*q)/r
end do
end do
end do
end subroutine potential3D
end module physics
Back to Main Page |