##// END OF EJS Templates
schain
RSS
fix line 461 y 665
fix line 461 y 665
rflores - - Load All Authors

File last commit:

r1647:4ac1ad84b857
r1791:e4984c469c9f isr
Show More
fitacf.f
621 lines | 14.3 KiB | text/x-fortran | FortranFixedLexer
/ schainf / Ffiles / fitacf.f
History | Source | Raw |Copy content |Copy permalink
rflores
add Ffiles
 r1601 subroutine guess(acf,tau,npts,zero,amin,te,tr)
c
c find zero crossing (zero), depth of minimum (amin), height of maximum
c
real acf(npts),tau(npts)
zero=0.0
amin=1.0
tmin=0.0
jmin=0
do i=npts,2,-1
if(acf(i)*acf(i-1).lt.0.0) then
zero=(tau(i-1)*acf(i)-tau(i)*acf(i-1))/(acf(i)-acf(i-1))
end if
if(acf(i).lt.amin) then
amin=acf(i)
jmin=i
end if
end do
if(jmin.gt.0) then
call parab1(tau(jmin-1),acf(jmin-1),a,b,c)
tmin=-b/(2.0*a)
amin=c+tmin*(b+tmin*a)
end if
tr=cdtr1(-amin)
te=czte1(zero*1000.0,tr)
return
end
subroutine parab1(x,y,a,b,c)
C-----
dimension x(3),y(3)
delta=x(1)-x(2)
a=(y(1)-2.*y(2)+y(3))/(2.*delta*delta)
b=(y(1)-y(2))/delta - a*(x(1)+x(2))
c=y(1)-a*x(1)*x(1)-b*x(1)
return
end
real function cdtr1(depth)
C-----convert depth to te/ti ratio
dimension tr(4)
C according to the curve published in farley et al 1967
c modified for 2004 conditions on axis
c data tr/7.31081,3.53286,5.92271,.174/
data tr/9.5,4.0,8.5,.3/
data nt/4/
cdtr1=tr(1)
do i=2,nt
cdtr1=cdtr1*depth + tr(i)
end do
return
end
real function czte1(zlag,tr)
C-----convert zero crossing point to te
C according to the curve published in farley et al 1967
c modified for 2004 conditions on axis
dimension dt(4)
c data dt/0.00945025,-0.0774338,.203626,.812397/,nd/4/
data dt/0.00945025,-0.0774338,.2,0.9/,nd/4/
data t0/1000./
tr1=min(abs(tr),5.)
if(zlag .eq. 0)then
czte1=1000000.
else
dt0=dt(1)
do i=2,nd
dt0=dt0*tr1 + dt(i)
end do
czte1=t0*(dt0/zlag)**2
end if
return
end
subroutine fit(wl,taup,rhop,covar,cinv,sigma2p,paramp,ebp,
& bfldp,alphap,densp,alt,time,nl,ifitp,ist)
c
c subroutine to fit measured ACF
c wavelength wl (m),lags taup (s), normalized acf rhop,
c experimental variances sigma2p, covariances covar,
c fit parameters params, error bars ebp, magnetic field bfldp (gauss),
c alphap B-field angle (radians), density densp (gcs),
c altitude alt (km), time (LT hours), nl lags
c ifitp determines which parameters are fit (see below)
c
real tol,pi,wl
parameter(nlmax=100,npmax=10,lwa=2000,tol=1.0e-5)
external fcn
real covar(nl,nl),cinv(nl,nl)
real ev(nlmax*nlmax),ap(nlmax*(nlmax+1)/2)
real fv1(nlmax),fv2(nlmax),det(2),w(nlmax)
real taup(nl),rhop(nl),sigma2p(nl),paramp(npmax)
real p2(npmax),ebp(npmax),alphap,bfldp,densp
real tau(nlmax),rho(nlmax),sigma2(nlmax),params(npmax)
real fvec(nlmax),wa(lwa),wa2(lwa),eb(npmax)
integer iwa(npmax),ifitp(npmax),ifit(npmax)
include 'fitter.h'
common /mode/imode
common/fitter/tau,rho,sigma2,params,ifit
common/trans/ev
C-----
c set ifit(1-5) to unity to fit:
c 1 normalization
c 2 Te
c 3 Ti
c 4 H+
c 5 He+
imode=2
c
pi=4.0*atan(1.0)
ak=2.0*pi/wl
wi(1)=16
wi(2)=1
wi(3)=4
nion=3
c
c invert covariances and find eigenvalues and eigenvectors
c
l = 0
do 20 j = 1, nl
do 10 i = 1, j
l = l + 1
ap(l) = covar(i,j)
10 continue
20 continue
call sppfa(ap,nl,info)
call sppdi(ap,nl,det,01)
l = 0
do 40 j = 1, nl
do 30 i = 1, j
l = l + 1
cinv(i,j)=ap(l)
cinv(j,i)=ap(l)
30 continue
40 continue
c
c transformation matrix is inverse (transpose) of eigenvectors
c
matz=1
call rs(nl,nl,cinv,w,matz,ev,fv1,fv2,ierr)
do i=1,nl
sigma2(i)=1.0/w(i)
end do
c
c extract desired parameters
c
iparm=0
do i=1,5
eb(i)=0.0
if(ifitp(i).eq.1) then
iparm=iparm+1
p2(iparm)=paramp(i)
end if
ifit(i)=ifitp(i)
params(i)=paramp(i)
end do
np=iparm
alpha=alphap
dens=densp
do i=1,nl
tau(i)=taup(i)
rho(i)=rhop(i)
sigma2p(i)=sigma2(i)
end do
bfld=bfldp
c
c no. equations is no. lags - do nlls fit
c
call lmdif1(fcn,nl,np,p2,fvec,tol,info,iwa,wa,lwa)
ist=info
c
c generate error bars here
c
call fdjac2(fcn,nl,np,p2,fvec,wa,nl,iflag,0.0e0,wa2)
do i=1,np
err=0.0
do j=1,nl
err=err+(wa(j+(i-1)*nl))**2
end do
if(err.gt.0.0) eb(i)=sqrt(err**-1)
end do
c reorder results
iparm=0
do i=1,5
if(ifit(i).eq.1) then
iparm=iparm+1
paramp(i)=p2(iparm)
ebp(i)=eb(iparm)
else
ebp(i)=0.0
paramp(i)=params(i)
end if
end do
9 continue
c write(*,*) dens
c write(*,*) te
return
end
subroutine fcn(m,n,x,fvec,iflag)
c
c provides m functions (fvec) in n variables (x) to minimize
c in least-squares sense
c
parameter(nlmax=100,npmax=10)
integer m,n,iflag
real fvec(m)
integer ifit(npmax)
real x(n),fv(nlmax),anorm,params(npmax),ev(nlmax*nlmax)
real tau(nlmax),rho(nlmax),sigma2(nlmax)
real chisq,acf
include 'fitter.h'
common/fitter/tau,rho,sigma2,params,ifit
common /errs/ chisq
common /trans/ev
c 1 0 0 0 0 fit zero lag (otherwise set to default)
c 0 1 0 0 0 fit Te (otherwise default)
c 0 0 1 0 0 fit Ti (otherwise set to Te)
c 0 0 0 1 0 fit H+ (otherwize default)
c 0 0 0 0 1 fit He+ (otherwise default)
8 continue
c
c ignore collisions
c
c write(*,*) "starting fcn"
ven=0.0
do i=1,3
vin(i)=0.0
end do
iparm=0
c
c zero lag normalization constant
c
if(ifit(1).eq.1) then
iparm=iparm+1
c=x(iparm)
else
c=params(1)
end if
c
c Te
c
if(ifit(2).eq.1) then
iparm=iparm+1
te=x(iparm)
else
te=params(2)
end if
c
c Ti - default is Te rather than initial value
c
if(ifit(3).eq.1) then
iparm=iparm+1
do i=1,3
ti(i)=x(iparm)
end do
else
do i=1,3
ti(i)=te
end do
params(3)=te
end if
c three-ion plasma
nion=3
c
c composition H+ first
c
fi(1)=1.0
if(ifit(4).eq.1) then
iparm=iparm+1
fi(2)=(x(iparm))
fi(1)=fi(1)-(x(iparm))
else
fi(2)=params(4)
fi(1)=fi(1)-params(4)
end if
c
c He+
c
if(ifit(5).eq.1) then
iparm=iparm+1
fi(3)=(x(iparm))
fi(1)=fi(1)-(x(iparm))
else
fi(3)=params(5)
fi(1)=fi(1)-params(5)
end if
c write(*,*) x,bfld,alpha,dens,ak,m
call gaussq(0.0,anorm)
anorm=anorm/c
if(anorm.eq.0.0.or.m.eq.0)return
chisq=0.0
do i=1,m
call gaussq(tau(i),acf)
fv(i)=(acf/anorm-rho(i))
end do
c
c transform to space where s2 is diagonal
c (are i and j transposed below?)
c
do i=1,m
fvec(i)=0.0
do j=1,m
fvec(i)=fvec(i)+fv(j)*ev(j+(i-1)*m)/sqrt(sigma2(i))
end do
chisq=chisq+fvec(i)**2
end do
chisq=chisq/float(m)
c write(*,*) fvec
c write(*,*) chisq
c stop
return
end
c
c not really fond of above code
c
complex function cj_ion(theta,psi)
c
real theta,phi,psi,alpha
complex z
complex*16 zz
z=zz(dcmplx(-theta,psi))
cj_ion=z*cmplx(0.0,-1.0)
return
end
complex function cj_electron(theta,phi,psi,alpha)
c
c Theta, phi, and psi are the normalized frequency, gyrofrequency,
c and collision frequency. Alpha is angle between wavevector and
c magnetic field in radians.
c
parameter(nterms=10)
real theta,phi,psi,alpha
real arg
complex cj,cy(0:nterms)
complex z
complex*16 zz
integer m,nz,ierr
integer imode
common/mode/imode
c
if(imode.eq.3) then
arg=0.5*(sin(alpha)/phi)**2
c calculate modified Bessel functions using amos library
call cbesi(cmplx(arg,0.0),0.0,2,nterms+1,cy,nz,ierr)
cj=cmplx(0.0,0.0)
do m=-nterms,nterms
z=zz(dcmplx(-(theta-float(m)*phi)/cos(alpha),
& psi/cos(alpha)))
cj=cj+z*cy(iabs(m))
end do
cj_electron=cj*cmplx(0.0,-1.0/cos(alpha))
else
z=zz(dcmplx(-theta/cos(alpha),psi/cos(alpha)))
cj_electron=z*cmplx(0.0,-1.0/cos(alpha))
cj_electron=cj_electron*(1.0-(sin(alpha)/phi)**2/2.0)
end if
return
end
complex function y_ion(theta,psi)
c
real theta,phi,psi,alpha
complex cj, cj_ion
cj=cj_ion(theta,psi)
y_ion=cj*cmplx(theta,-psi)+cmplx(0.0,1.0)
y_ion=y_ion/(1-psi*cj)
return
end
complex function y_electron(theta,phi,psi,alpha)
c
real theta,phi,psi,alpha
complex cj, cj_electron
cj=cj_electron(theta,phi,psi,alpha)
y_electron=cj*cmplx(theta,-psi)+cmplx(0.0,1.0)
y_electron=y_electron/(1.0-psi*cj)
return
end
real function spect1(omega)
c
c Function to generate IS ion-line spectrum
c Provisions for nimax ions
c No provisions for drifts
c
c imode=1 Farley B-field treatment for electrons
c 2 use Mike Sulzer's model for ye
c 3 calculate ye using sums of Bessel functions
c
include 'fitter.h'
real omega,thetae,thetai(nimax),psie,psii(nimax),phi,p
real tr(nimax),vti(nimax),vte,omegae
real bk,em,e,dlf,dl,pi
real alpha2,densmks,freq
complex ye,yed,yet,yi(nimax),sum1,sum2,sumdl
complex y_ion, y_electron, y_esum
integer i,j,k,imode
common/mode/ imode
data bk,em,e,dlf/1.38e-23,9.1e-31,1.6e-19,4772.9/
c fi is ion fraction, wi is ion atomic weight
c dens is electron density (cgs), alpha is angle
c between k and the magnetic field (radians)
pi=4.0*atan(1.0)
if(omega.eq.0.0) omega=1.0
omegae=e*bfld*1.0e-4/em
densmks=dens*1.0e6
sum1=0.0
sum2=0.0
do i=1,nion
tr(i)=te/ti(i)
vti(i)=sqrt(bk*ti(i)/(1.67e-27*wi(i)))
thetai(i)=(omega/ak)/(sqrt(2.0)*vti(i))
psii(i)=(vin(i)/ak)/(sqrt(2.0)*vti(i))
yi(i)=y_ion(thetai(i),psii(i))
sum1=sum1+fi(i)*tr(i)*yi(i)
sum2=sum2+fi(i)*yi(i)
end do
dl=ak**2*dlf*te/densmks
rflores
test 5
 r1640 c write(*,fmt='("Before YE")')
c write(*,*) imode
c call exit
rflores
add Ffiles
 r1601
if(imode.eq.1.or.imode.eq.3) then
vte=sqrt(bk*te/em)
thetae=(omega/ak)/(sqrt(2.0)*vte)
phi=(omegae/ak)/(sqrt(2.0)*vte)
psie=(ven/ak)/(sqrt(2.0)*vte)
ye=y_electron(thetae,phi,psie,alpha)
rflores
test 5
 r1640 write(*,fmt='("AFTER YE")')
c call exit
rflores
add Ffiles
 r1601
else if(imode.eq.2) then
c
c use Mike Sulzer's library here: alpha2 is angle off perp in degrees
c
freq=omega/(2.0*pi)
alpha2=abs(pi/2.0-alpha)*180.0/pi
c write(*,*) "ye: ", ye
call collision(densmks, te, freq, alpha2, ye)
rflores
test 5
 r1640 c write(*,fmt='("AFTER COLLISION")')
c call exit
rflores
add Ffiles
 r1601 ye=ye*omega+cmplx(0.0,1.0)
end if
yed=ye+cmplx(0.0,dl)
p=(cabs(ye))**2*real(sum2)+cabs(sum1+cmplx(0.0,dl))**2*real(ye)
p=p/(cabs(yed+sum1))**2
spect1=p*2.0e0/(omega*pi)
rflores
test 5
 r1640 write(*,*) "spect1:",spect1
rflores
add Ffiles
 r1601 return
end
rflores
fixed fitacf.f
 r1647 subroutine acf2(wl, tau, te1, ti1, fi1, ven1, vin1, wi1, nion1,
rflores
bug fixes
 r1645 & alpha1, dens1, bfld1, acf)
rflores
add Ffiles
 r1601 c
c computes autocorrelation function for given plasma parameters
c by integrating real spectrum
c tau in sec., alpha in radians, density in cgs, bfield in cgs
c scattering wavelength (wl) in meters
c
include 'fitter.h'
real wl,tau,te1,ti1(nion1),fi1(nion1),ven1,vin1(nion1),alpha1,
& dens1,bfld1,acf
real pi
integer nion1
integer wi1(nion1)
rflores
test 5
 r1640 integer i,j,k,imode
common /mode/imode
rflores
add Ffiles
 r1601 c
rflores
test 5
 r1640 c write(*,*) "INITIAL acf:",wl,tau,te1,ti1,fi1,ven1,vin1,wi1,alpha1
write(*,*) "INITIAL acf:",dens1, bfld1, acf, nion1
c write(*,fmt='("INIT")')
rflores
add Ffiles
 r1601 pi=4.0*atan(1.0)
c
c copy arguments to common block
c
ak=2.0*pi/wl
imode=2
rflores
test 5
 r1640 write(*,*) "imode:",imode
rflores
add Ffiles
 r1601
nion=nion1
alpha=alpha1
te=te1
ven=ven1
rflores
test 5
 r1640 c write(*,fmt='("INIT2")')
rflores
add Ffiles
 r1601 do i=1,nion
rflores
test 5
 r1640 c write(*,fmt='("INIT2.5")')
rflores
add Ffiles
 r1601 ti(i)=ti1(i)
fi(i)=fi1(i)
vin(i)=vin1(i)
wi(i)=wi1(i)
end do
rflores
test 5
 r1640 c write(*,fmt='("INIT3")')
rflores
add Ffiles
 r1601 dens=dens1
bfld=bfld1
c write(*,*) wl,alpha1,bfld1,dens
c call exit
rflores
test 5
 r1640 c write(*,fmt='("Before Gauss")')
rflores
fitacf.f update
 r1636 call gaussq(tau,acf)
rflores
add Ffiles
 r1601
rflores
test 5
 r1640 write(*,*) "FINAL acf:",acf
c write(*,fmt='("After Gauss")')
rflores
add Ffiles
 r1601 return
end
subroutine gaussq(tau,acf)
c
c Computes cosine or sine transform of given real function
c Uses quadpack, tau in sec.
c
real a,abserr,epsabs,acf,tau,work,chebmo
integer ier,integr,iwork,last,leniw,lenw,limit,limlst,
* lst,maxp1,neval
integer knorm,ksave,momcom,nrmom
dimension iwork(300),work(1625),chebmo(61,25)
external spect1
include 'fitter.h'
c
c write(*,*) "inside gaussq"
a = 0.0
b = 2.5e4*ak ! upper integration limit open to debate (was 1.5, now 2.5)
integr = 1
epsabs = 1.0e-4
limlst = 50
limit = 100
leniw = limit*2+limlst
maxp1 = 61
lenw = leniw*2+maxp1*25
c write(*,*) leniw,lenw
nrmom=0
ksave=0
momcom=0
rflores
test 5
 r1640 write(*,*) "Before qc25f:",acf,imode
rflores
add Ffiles
 r1601 c much faster, more robust
c write(*,*) "acf_in: ",acf
call qc25f(spect1,a,b,tau,integr,nrmom,maxp1,ksave,acf,
& abserr,neval,resabs,resasc,momcom,chebmo)
c write(*,*) "acf_out: ",acf
c call qawf(spect1,a,tau,integr,epsabs,acf,abserr,neval,
c & ier,limlst,lst,leniw,maxp1,lenw,iwork,work)
rflores
test 5
 r1640 write(*,*) "After qc25f:",acf
c call exit
rflores
add Ffiles
 r1601 return
end