program dielec
C COMPUTE DIELECTRIC FROM TIME SERIES OF BOX DIPOLE MAGNITUDE
C use Coulombs, meters, Joules
integer i,j,k, nt,mxt
character*80 fn
parameter (mxt=10000)
real*8 dp(mxt), d2(mxt), ti(mxt),ed(mxt), bx(mxt),by(mxt),bz(mxt)
real*8 adp,ad2, pi,tem, tv,av, mx,my,mz, rdi,rd2, sd2
C KB == Boltzmann, DC == Debye to C.m, EPS0 == vac permittivity
REAL*8 KB, DC, EPS0
PARAMETER (KB=1.3805d-23, DC=3.336d-30, EPS0=8.85419d-12)
pi = acos(-1.d0)
tem = 293.15d0
k=iargc()
if (k.ne.2) then
write(6,*) 'dielec.ix DIPOLE.DAT CUMDIEL.DAT'
write(6,*) ' DIPOLE.DAT input; time x y z edge dipole_mag'
write(6,*) ' CUMDIEL.DAT output; time cum_dielec'
stop 'and try again'
else
call getarg(1,fn)
open(1,file=fn,form='formatted')
call getarg(2,fn)
open(2,file=fn,form='formatted')
endif
nt = 1
10 read(1,*,end=20,err=20) ti(nt),bx(nt),by(nt),bz(nt),ed(nt),dp(nt)
nt = nt + 1
if (nt.le.mxt) goto 10
20 continue
if (nt.lt.mxt) nt = nt-1
fac = (dc*dc) / (3.d0*kb*tem*eps0)
write (6,*) 'Scale factor = ', fac
tv = 0.d0
mx = 0.d0
my = 0.d0
mz = 0.d0
sd2 = 0.d0
do k=1,nt
mx = mx + bx(k)
my = my + by(k)
mz = mz + bz(k)
rdi = sqrt(mx**2 + my**2 + mz**2)
sd2 = sd2 + bx(k)**2 + by(k)**2 + bz(k)**2
tv = tv + (1.d-10*ed(k))**3
av = tv / dble(k)
adp = rdi / dble(k)
ad2 = sd2 / dble(k)
write (2,'(2F11.1)') ti(k),1.d0 + fac*(ad2-adp**2)/av
enddo
write (6,*) 'Final eps=', 1.d0 + fac*(ad2-adp**2)/av,
1 ' < M**2 > =',ad2, ' < M > =',adp
stop
end
