schain Commit - r1636:43055b8aca7f · Jicamarca Repository

1

subroutine guess(acf,tau,npts,zero,amin,te,tr)

1

subroutine guess(acf,tau,npts,zero,amin,te,tr)

2

c

2

c

3

c find zero crossing (zero), depth of minimum (amin), height of maximum

3

c find zero crossing (zero), depth of minimum (amin), height of maximum

4

c

4

c

5

real acf(npts),tau(npts)

5

real acf(npts),tau(npts)

6

7

zero=0.0

7

zero=0.0

8

amin=1.0

8

amin=1.0

9

tmin=0.0

9

tmin=0.0

10

jmin=0

10

jmin=0

11

12

do i=npts,2,-1

12

do i=npts,2,-1

13

if(acf(i)*acf(i-1).lt.0.0) then

13

if(acf(i)*acf(i-1).lt.0.0) then

14

zero=(tau(i-1)*acf(i)-tau(i)*acf(i-1))/(acf(i)-acf(i-1))

14

zero=(tau(i-1)*acf(i)-tau(i)*acf(i-1))/(acf(i)-acf(i-1))

15

end if

15

end if

16

if(acf(i).lt.amin) then

16

if(acf(i).lt.amin) then

17

amin=acf(i)

17

amin=acf(i)

18

jmin=i

18

jmin=i

19

end if

19

end if

20

end do

20

end do

21

22

if(jmin.gt.0) then

22

if(jmin.gt.0) then

23

call parab1(tau(jmin-1),acf(jmin-1),a,b,c)

23

call parab1(tau(jmin-1),acf(jmin-1),a,b,c)

24

tmin=-b/(2.0*a)

24

tmin=-b/(2.0*a)

25

amin=c+tmin*(b+tmin*a)

25

amin=c+tmin*(b+tmin*a)

26

end if

26

end if

27

28

tr=cdtr1(-amin)

28

tr=cdtr1(-amin)

29

te=czte1(zero*1000.0,tr)

29

te=czte1(zero*1000.0,tr)

30

return

30

return

31

end

31

end

32

33

subroutine parab1(x,y,a,b,c)

33

subroutine parab1(x,y,a,b,c)

34

C-----

34

C-----

35

dimension x(3),y(3)

35

dimension x(3),y(3)

36

delta=x(1)-x(2)

36

delta=x(1)-x(2)

37

a=(y(1)-2.*y(2)+y(3))/(2.*delta*delta)

37

a=(y(1)-2.*y(2)+y(3))/(2.*delta*delta)

38

b=(y(1)-y(2))/delta - a*(x(1)+x(2))

38

b=(y(1)-y(2))/delta - a*(x(1)+x(2))

39

c=y(1)-a*x(1)*x(1)-b*x(1)

39

c=y(1)-a*x(1)*x(1)-b*x(1)

40

return

40

return

41

end

41

end

42

43

real function cdtr1(depth)

43

real function cdtr1(depth)

44

C-----convert depth to te/ti ratio

44

C-----convert depth to te/ti ratio

45

dimension tr(4)

45

dimension tr(4)

46

C according to the curve published in farley et al 1967

46

C according to the curve published in farley et al 1967

47

c modified for 2004 conditions on axis

47

c modified for 2004 conditions on axis

48

c data tr/7.31081,3.53286,5.92271,.174/

48

c data tr/7.31081,3.53286,5.92271,.174/

49

data tr/9.5,4.0,8.5,.3/

49

data tr/9.5,4.0,8.5,.3/

50

data nt/4/

50

data nt/4/

51

cdtr1=tr(1)

51

cdtr1=tr(1)

52

do i=2,nt

52

do i=2,nt

53

cdtr1=cdtr1*depth + tr(i)

53

cdtr1=cdtr1*depth + tr(i)

54

end do

54

end do

55

return

55

return

56

end

56

end

57

58

real function czte1(zlag,tr)

58

real function czte1(zlag,tr)

59

C-----convert zero crossing point to te

59

C-----convert zero crossing point to te

60

C according to the curve published in farley et al 1967

60

C according to the curve published in farley et al 1967

61

c modified for 2004 conditions on axis

61

c modified for 2004 conditions on axis

62

dimension dt(4)

62

dimension dt(4)

63

c data dt/0.00945025,-0.0774338,.203626,.812397/,nd/4/

63

c data dt/0.00945025,-0.0774338,.203626,.812397/,nd/4/

64

data dt/0.00945025,-0.0774338,.2,0.9/,nd/4/

64

data dt/0.00945025,-0.0774338,.2,0.9/,nd/4/

65

data t0/1000./

65

data t0/1000./

66

tr1=min(abs(tr),5.)

66

tr1=min(abs(tr),5.)

67

if(zlag .eq. 0)then

67

if(zlag .eq. 0)then

68

czte1=1000000.

68

czte1=1000000.

69

else

69

else

70

dt0=dt(1)

70

dt0=dt(1)

71

do i=2,nd

71

do i=2,nd

72

dt0=dt0*tr1 + dt(i)

72

dt0=dt0*tr1 + dt(i)

73

end do

73

end do

74

czte1=t0*(dt0/zlag)**2

74

czte1=t0*(dt0/zlag)**2

75

end if

75

end if

76

return

76

return

77

end

77

end

78

79

subroutine fit(wl,taup,rhop,covar,cinv,sigma2p,paramp,ebp,

79

subroutine fit(wl,taup,rhop,covar,cinv,sigma2p,paramp,ebp,

80

& bfldp,alphap,densp,alt,time,nl,ifitp,ist)

80

& bfldp,alphap,densp,alt,time,nl,ifitp,ist)

81

c

81

c

82

c subroutine to fit measured ACF

82

c subroutine to fit measured ACF

83

c wavelength wl (m),lags taup (s), normalized acf rhop,

83

c wavelength wl (m),lags taup (s), normalized acf rhop,

84

c experimental variances sigma2p, covariances covar,

84

c experimental variances sigma2p, covariances covar,

85

c fit parameters params, error bars ebp, magnetic field bfldp (gauss),

85

c fit parameters params, error bars ebp, magnetic field bfldp (gauss),

86

c alphap B-field angle (radians), density densp (gcs),

86

c alphap B-field angle (radians), density densp (gcs),

87

c altitude alt (km), time (LT hours), nl lags

87

c altitude alt (km), time (LT hours), nl lags

88

c ifitp determines which parameters are fit (see below)

88

c ifitp determines which parameters are fit (see below)

89

c

89

c

90

real tol,pi,wl

90

real tol,pi,wl

91

parameter(nlmax=100,npmax=10,lwa=2000,tol=1.0e-5)

91

parameter(nlmax=100,npmax=10,lwa=2000,tol=1.0e-5)

92

external fcn

92

external fcn

93

94

real covar(nl,nl),cinv(nl,nl)

94

real covar(nl,nl),cinv(nl,nl)

95

real ev(nlmax*nlmax),ap(nlmax*(nlmax+1)/2)

95

real ev(nlmax*nlmax),ap(nlmax*(nlmax+1)/2)

96

real fv1(nlmax),fv2(nlmax),det(2),w(nlmax)

96

real fv1(nlmax),fv2(nlmax),det(2),w(nlmax)

97

98

real taup(nl),rhop(nl),sigma2p(nl),paramp(npmax)

98

real taup(nl),rhop(nl),sigma2p(nl),paramp(npmax)

99

real p2(npmax),ebp(npmax),alphap,bfldp,densp

99

real p2(npmax),ebp(npmax),alphap,bfldp,densp

100

real tau(nlmax),rho(nlmax),sigma2(nlmax),params(npmax)

100

real tau(nlmax),rho(nlmax),sigma2(nlmax),params(npmax)

101

real fvec(nlmax),wa(lwa),wa2(lwa),eb(npmax)

101

real fvec(nlmax),wa(lwa),wa2(lwa),eb(npmax)

102

103

integer iwa(npmax),ifitp(npmax),ifit(npmax)

103

integer iwa(npmax),ifitp(npmax),ifit(npmax)

104

105

include 'fitter.h'

105

include 'fitter.h'

106

common /mode/imode

106

common /mode/imode

107

108

common/fitter/tau,rho,sigma2,params,ifit

108

common/fitter/tau,rho,sigma2,params,ifit

109

common/trans/ev

109

common/trans/ev

110

C-----

110

C-----

111

112

c set ifit(1-5) to unity to fit:

112

c set ifit(1-5) to unity to fit:

113

c 1 normalization

113

c 1 normalization

114

c 2 Te

114

c 2 Te

115

c 3 Ti

115

c 3 Ti

116

c 4 H+

116

c 4 H+

117

c 5 He+

117

c 5 He+

118

119

imode=2

119

imode=2

120

c

120

c

121

pi=4.0*atan(1.0)

121

pi=4.0*atan(1.0)

122

ak=2.0*pi/wl

122

ak=2.0*pi/wl

123

124

wi(1)=16

124

wi(1)=16

125

wi(2)=1

125

wi(2)=1

126

wi(3)=4

126

wi(3)=4

127

nion=3

127

nion=3

128

129

c

129

c

130

c invert covariances and find eigenvalues and eigenvectors

130

c invert covariances and find eigenvalues and eigenvectors

131

c

131

c

132

l = 0

132

l = 0

133

do 20 j = 1, nl

133

do 20 j = 1, nl

134

do 10 i = 1, j

134

do 10 i = 1, j

135

l = l + 1

135

l = l + 1

136

ap(l) = covar(i,j)

136

ap(l) = covar(i,j)

137

10 continue

137

10 continue

138

20 continue

138

20 continue

139

140

call sppfa(ap,nl,info)

140

call sppfa(ap,nl,info)

141

call sppdi(ap,nl,det,01)

141

call sppdi(ap,nl,det,01)

142

143

l = 0

143

l = 0

144

do 40 j = 1, nl

144

do 40 j = 1, nl

145

do 30 i = 1, j

145

do 30 i = 1, j

146

l = l + 1

146

l = l + 1

147

cinv(i,j)=ap(l)

147

cinv(i,j)=ap(l)

148

cinv(j,i)=ap(l)

148

cinv(j,i)=ap(l)

149

30 continue

149

30 continue

150

40 continue

150

40 continue

151

152

c

152

c

153

c transformation matrix is inverse (transpose) of eigenvectors

153

c transformation matrix is inverse (transpose) of eigenvectors

154

c

154

c

155

156

matz=1

156

matz=1

157

call rs(nl,nl,cinv,w,matz,ev,fv1,fv2,ierr)

157

call rs(nl,nl,cinv,w,matz,ev,fv1,fv2,ierr)

158

159

do i=1,nl

159

do i=1,nl

160

sigma2(i)=1.0/w(i)

160

sigma2(i)=1.0/w(i)

161

end do

161

end do

162

163

c

163

c

164

c extract desired parameters

164

c extract desired parameters

165

c

165

c

166

iparm=0

166

iparm=0

167

do i=1,5

167

do i=1,5

168

eb(i)=0.0

168

eb(i)=0.0

169

if(ifitp(i).eq.1) then

169

if(ifitp(i).eq.1) then

170

iparm=iparm+1

170

iparm=iparm+1

171

p2(iparm)=paramp(i)

171

p2(iparm)=paramp(i)

172

end if

172

end if

173

ifit(i)=ifitp(i)

173

ifit(i)=ifitp(i)

174

params(i)=paramp(i)

174

params(i)=paramp(i)

175

end do

175

end do

176

np=iparm

176

np=iparm

177

178

alpha=alphap

178

alpha=alphap

179

dens=densp

179

dens=densp

180

do i=1,nl

180

do i=1,nl

181

tau(i)=taup(i)

181

tau(i)=taup(i)

182

rho(i)=rhop(i)

182

rho(i)=rhop(i)

183

sigma2p(i)=sigma2(i)

183

sigma2p(i)=sigma2(i)

184

end do

184

end do

185

186

bfld=bfldp

186

bfld=bfldp

187

188

c

188

c

189

c no. equations is no. lags - do nlls fit

189

c no. equations is no. lags - do nlls fit

190

c

190

c

191

192

call lmdif1(fcn,nl,np,p2,fvec,tol,info,iwa,wa,lwa)

192

call lmdif1(fcn,nl,np,p2,fvec,tol,info,iwa,wa,lwa)

193

ist=info

193

ist=info

194

195

c

195

c

196

c generate error bars here

196

c generate error bars here

197

c

197

c

198

call fdjac2(fcn,nl,np,p2,fvec,wa,nl,iflag,0.0e0,wa2)

198

call fdjac2(fcn,nl,np,p2,fvec,wa,nl,iflag,0.0e0,wa2)

199

200

do i=1,np

200

do i=1,np

201

err=0.0

201

err=0.0

202

do j=1,nl

202

do j=1,nl

203

err=err+(wa(j+(i-1)*nl))**2

203

err=err+(wa(j+(i-1)*nl))**2

204

end do

204

end do

205

if(err.gt.0.0) eb(i)=sqrt(err**-1)

205

if(err.gt.0.0) eb(i)=sqrt(err**-1)

206

end do

206

end do

207

208

c reorder results

208

c reorder results

209

iparm=0

209

iparm=0

210

do i=1,5

210

do i=1,5

211

if(ifit(i).eq.1) then

211

if(ifit(i).eq.1) then

212

iparm=iparm+1

212

iparm=iparm+1

213

paramp(i)=p2(iparm)

213

paramp(i)=p2(iparm)

214

ebp(i)=eb(iparm)

214

ebp(i)=eb(iparm)

215

else

215

else

216

ebp(i)=0.0

216

ebp(i)=0.0

217

paramp(i)=params(i)

217

paramp(i)=params(i)

218

end if

218

end if

219

end do

219

end do

220

221

9 continue

221

9 continue

222

c write(*,*) dens

222

c write(*,*) dens

223

c write(*,*) te

223

c write(*,*) te

224

225

return

225

return

226

end

226

end

227

228

229

subroutine fcn(m,n,x,fvec,iflag)

229

subroutine fcn(m,n,x,fvec,iflag)

230

c

230

c

231

c provides m functions (fvec) in n variables (x) to minimize

231

c provides m functions (fvec) in n variables (x) to minimize

232

c in least-squares sense

232

c in least-squares sense

233

c

233

c

234

235

parameter(nlmax=100,npmax=10)

235

parameter(nlmax=100,npmax=10)

236

integer m,n,iflag

236

integer m,n,iflag

237

real fvec(m)

237

real fvec(m)

238

239

integer ifit(npmax)

239

integer ifit(npmax)

240

real x(n),fv(nlmax),anorm,params(npmax),ev(nlmax*nlmax)

240

real x(n),fv(nlmax),anorm,params(npmax),ev(nlmax*nlmax)

241

real tau(nlmax),rho(nlmax),sigma2(nlmax)

241

real tau(nlmax),rho(nlmax),sigma2(nlmax)

242

real chisq,acf

242

real chisq,acf

243

244

include 'fitter.h'

244

include 'fitter.h'

245

246

common/fitter/tau,rho,sigma2,params,ifit

246

common/fitter/tau,rho,sigma2,params,ifit

247

common /errs/ chisq

247

common /errs/ chisq

248

common /trans/ev

248

common /trans/ev

249

250

c 1 0 0 0 0 fit zero lag (otherwise set to default)

250

c 1 0 0 0 0 fit zero lag (otherwise set to default)

251

c 0 1 0 0 0 fit Te (otherwise default)

251

c 0 1 0 0 0 fit Te (otherwise default)

252

c 0 0 1 0 0 fit Ti (otherwise set to Te)

252

c 0 0 1 0 0 fit Ti (otherwise set to Te)

253

c 0 0 0 1 0 fit H+ (otherwize default)

253

c 0 0 0 1 0 fit H+ (otherwize default)

254

c 0 0 0 0 1 fit He+ (otherwise default)

254

c 0 0 0 0 1 fit He+ (otherwise default)

255

256

8 continue

256

8 continue

257

258

c

258

c

259

c ignore collisions

259

c ignore collisions

260

c

260

c

261

c write(*,*) "starting fcn"

261

c write(*,*) "starting fcn"

262

ven=0.0

262

ven=0.0

263

do i=1,3

263

do i=1,3

264

vin(i)=0.0

264

vin(i)=0.0

265

end do

265

end do

266

267

iparm=0

267

iparm=0

268

269

c

269

c

270

c zero lag normalization constant

270

c zero lag normalization constant

271

c

271

c

272

if(ifit(1).eq.1) then

272

if(ifit(1).eq.1) then

273

iparm=iparm+1

273

iparm=iparm+1

274

c=x(iparm)

274

c=x(iparm)

275

else

275

else

276

c=params(1)

276

c=params(1)

277

end if

277

end if

278

c

278

c

279

c Te

279

c Te

280

c

280

c

281

if(ifit(2).eq.1) then

281

if(ifit(2).eq.1) then

282

iparm=iparm+1

282

iparm=iparm+1

283

te=x(iparm)

283

te=x(iparm)

284

else

284

else

285

te=params(2)

285

te=params(2)

286

end if

286

end if

287

c

287

c

288

c Ti - default is Te rather than initial value

288

c Ti - default is Te rather than initial value

289

c

289

c

290

if(ifit(3).eq.1) then

290

if(ifit(3).eq.1) then

291

iparm=iparm+1

291

iparm=iparm+1

292

do i=1,3

292

do i=1,3

293

ti(i)=x(iparm)

293

ti(i)=x(iparm)

294

end do

294

end do

295

else

295

else

296

do i=1,3

296

do i=1,3

297

ti(i)=te

297

ti(i)=te

298

end do

298

end do

299

params(3)=te

299

params(3)=te

300

end if

300

end if

301

302

c three-ion plasma

302

c three-ion plasma

303

304

nion=3

304

nion=3

305

c

305

c

306

c composition H+ first

306

c composition H+ first

307

c

307

c

308

fi(1)=1.0

308

fi(1)=1.0

309

if(ifit(4).eq.1) then

309

if(ifit(4).eq.1) then

310

iparm=iparm+1

310

iparm=iparm+1

311

fi(2)=(x(iparm))

311

fi(2)=(x(iparm))

312

fi(1)=fi(1)-(x(iparm))

312

fi(1)=fi(1)-(x(iparm))

313

else

313

else

314

fi(2)=params(4)

314

fi(2)=params(4)

315

fi(1)=fi(1)-params(4)

315

fi(1)=fi(1)-params(4)

316

end if

316

end if

317

c

317

c

318

c He+

318

c He+

319

c

319

c

320

if(ifit(5).eq.1) then

320

if(ifit(5).eq.1) then

321

iparm=iparm+1

321

iparm=iparm+1

322

fi(3)=(x(iparm))

322

fi(3)=(x(iparm))

323

fi(1)=fi(1)-(x(iparm))

323

fi(1)=fi(1)-(x(iparm))

324

else

324

else

325

fi(3)=params(5)

325

fi(3)=params(5)

326

fi(1)=fi(1)-params(5)

326

fi(1)=fi(1)-params(5)

327

end if

327

end if

328

329

c write(*,*) x,bfld,alpha,dens,ak,m

329

c write(*,*) x,bfld,alpha,dens,ak,m

330

331

call gaussq(0.0,anorm)

331

call gaussq(0.0,anorm)

332

anorm=anorm/c

332

anorm=anorm/c

333

334

if(anorm.eq.0.0.or.m.eq.0)return

334

if(anorm.eq.0.0.or.m.eq.0)return

335

chisq=0.0

335

chisq=0.0

336

do i=1,m

336

do i=1,m

337

call gaussq(tau(i),acf)

337

call gaussq(tau(i),acf)

338

fv(i)=(acf/anorm-rho(i))

338

fv(i)=(acf/anorm-rho(i))

339

end do

339

end do

340

341

c

341

c

342

c transform to space where s2 is diagonal

342

c transform to space where s2 is diagonal

343

c (are i and j transposed below?)

343

c (are i and j transposed below?)

344

c

344

c

345

do i=1,m

345

do i=1,m

346

fvec(i)=0.0

346

fvec(i)=0.0

347

do j=1,m

347

do j=1,m

348

fvec(i)=fvec(i)+fv(j)*ev(j+(i-1)*m)/sqrt(sigma2(i))

348

fvec(i)=fvec(i)+fv(j)*ev(j+(i-1)*m)/sqrt(sigma2(i))

349

end do

349

end do

350

chisq=chisq+fvec(i)**2

350

chisq=chisq+fvec(i)**2

351

end do

351

end do

352

353

chisq=chisq/float(m)

353

chisq=chisq/float(m)

354

355

c write(*,*) fvec

355

c write(*,*) fvec

356

c write(*,*) chisq

356

c write(*,*) chisq

357

c stop

357

c stop

358

359

return

359

return

360

end

360

end

361

362

c

362

c

363

c not really fond of above code

363

c not really fond of above code

364

c

364

c

365

366

367

complex function cj_ion(theta,psi)

367

complex function cj_ion(theta,psi)

368

c

368

c

369

real theta,phi,psi,alpha

369

real theta,phi,psi,alpha

370

complex z

370

complex z

371

complex*16 zz

371

complex*16 zz

372

z=zz(dcmplx(-theta,psi))

372

z=zz(dcmplx(-theta,psi))

373

cj_ion=z*cmplx(0.0,-1.0)

373

cj_ion=z*cmplx(0.0,-1.0)

374

return

374

return

375

end

375

end

376

377

complex function cj_electron(theta,phi,psi,alpha)

377

complex function cj_electron(theta,phi,psi,alpha)

378

c

378

c

379

c Theta, phi, and psi are the normalized frequency, gyrofrequency,

379

c Theta, phi, and psi are the normalized frequency, gyrofrequency,

380

c and collision frequency. Alpha is angle between wavevector and

380

c and collision frequency. Alpha is angle between wavevector and

381

c magnetic field in radians.

381

c magnetic field in radians.

382

c

382

c

383

parameter(nterms=10)

383

parameter(nterms=10)

384

real theta,phi,psi,alpha

384

real theta,phi,psi,alpha

385

real arg

385

real arg

386

complex cj,cy(0:nterms)

386

complex cj,cy(0:nterms)

387

complex z

387

complex z

388

complex*16 zz

388

complex*16 zz

389

integer m,nz,ierr

389

integer m,nz,ierr

390

integer imode

390

integer imode

391

common/mode/imode

391

common/mode/imode

392

c

392

c

393

if(imode.eq.3) then

393

if(imode.eq.3) then

394

395

arg=0.5*(sin(alpha)/phi)**2

395

arg=0.5*(sin(alpha)/phi)**2

396

397

c calculate modified Bessel functions using amos library

397

c calculate modified Bessel functions using amos library

398

399

call cbesi(cmplx(arg,0.0),0.0,2,nterms+1,cy,nz,ierr)

399

call cbesi(cmplx(arg,0.0),0.0,2,nterms+1,cy,nz,ierr)

400

401

cj=cmplx(0.0,0.0)

401

cj=cmplx(0.0,0.0)

402

403

do m=-nterms,nterms

403

do m=-nterms,nterms

404

z=zz(dcmplx(-(theta-float(m)*phi)/cos(alpha),

404

z=zz(dcmplx(-(theta-float(m)*phi)/cos(alpha),

405

& psi/cos(alpha)))

405

& psi/cos(alpha)))

406

cj=cj+z*cy(iabs(m))

406

cj=cj+z*cy(iabs(m))

407

end do

407

end do

408

409

cj_electron=cj*cmplx(0.0,-1.0/cos(alpha))

409

cj_electron=cj*cmplx(0.0,-1.0/cos(alpha))

410

411

else

411

else

412

413

z=zz(dcmplx(-theta/cos(alpha),psi/cos(alpha)))

413

z=zz(dcmplx(-theta/cos(alpha),psi/cos(alpha)))

414

cj_electron=z*cmplx(0.0,-1.0/cos(alpha))

414

cj_electron=z*cmplx(0.0,-1.0/cos(alpha))

415

cj_electron=cj_electron*(1.0-(sin(alpha)/phi)**2/2.0)

415

cj_electron=cj_electron*(1.0-(sin(alpha)/phi)**2/2.0)

416

417

end if

417

end if

418

419

return

419

return

420

end

420

end

421

422

complex function y_ion(theta,psi)

422

complex function y_ion(theta,psi)

423

c

423

c

424

real theta,phi,psi,alpha

424

real theta,phi,psi,alpha

425

complex cj, cj_ion

425

complex cj, cj_ion

426

cj=cj_ion(theta,psi)

426

cj=cj_ion(theta,psi)

427

y_ion=cj*cmplx(theta,-psi)+cmplx(0.0,1.0)

427

y_ion=cj*cmplx(theta,-psi)+cmplx(0.0,1.0)

428

y_ion=y_ion/(1-psi*cj)

428

y_ion=y_ion/(1-psi*cj)

429

return

429

return

430

end

430

end

431

432

complex function y_electron(theta,phi,psi,alpha)

432

complex function y_electron(theta,phi,psi,alpha)

433

c

433

c

434

real theta,phi,psi,alpha

434

real theta,phi,psi,alpha

435

complex cj, cj_electron

435

complex cj, cj_electron

436

cj=cj_electron(theta,phi,psi,alpha)

436

cj=cj_electron(theta,phi,psi,alpha)

437

y_electron=cj*cmplx(theta,-psi)+cmplx(0.0,1.0)

437

y_electron=cj*cmplx(theta,-psi)+cmplx(0.0,1.0)

438

y_electron=y_electron/(1.0-psi*cj)

438

y_electron=y_electron/(1.0-psi*cj)

439

return

439

return

440

end

440

end

441

442

real function spect1(omega)

442

real function spect1(omega)

443

c

443

c

444

c Function to generate IS ion-line spectrum

444

c Function to generate IS ion-line spectrum

445

c Provisions for nimax ions

445

c Provisions for nimax ions

446

c No provisions for drifts

446

c No provisions for drifts

447

c

447

c

448

c imode=1 Farley B-field treatment for electrons

448

c imode=1 Farley B-field treatment for electrons

449

c 2 use Mike Sulzer's model for ye

449

c 2 use Mike Sulzer's model for ye

450

c 3 calculate ye using sums of Bessel functions

450

c 3 calculate ye using sums of Bessel functions

451

c

451

c

452

include 'fitter.h'

452

include 'fitter.h'

453

454

real omega,thetae,thetai(nimax),psie,psii(nimax),phi,p

454

real omega,thetae,thetai(nimax),psie,psii(nimax),phi,p

455

real tr(nimax),vti(nimax),vte,omegae

455

real tr(nimax),vti(nimax),vte,omegae

456

real bk,em,e,dlf,dl,pi

456

real bk,em,e,dlf,dl,pi

457

real alpha2,densmks,freq

457

real alpha2,densmks,freq

458

complex ye,yed,yet,yi(nimax),sum1,sum2,sumdl

458

complex ye,yed,yet,yi(nimax),sum1,sum2,sumdl

459

complex y_ion, y_electron, y_esum

459

complex y_ion, y_electron, y_esum

460

integer i,j,k,imode

460

integer i,j,k,imode

461

common/mode/ imode

461

common/mode/ imode

462

data bk,em,e,dlf/1.38e-23,9.1e-31,1.6e-19,4772.9/

462

data bk,em,e,dlf/1.38e-23,9.1e-31,1.6e-19,4772.9/

463

464

c fi is ion fraction, wi is ion atomic weight

464

c fi is ion fraction, wi is ion atomic weight

465

c dens is electron density (cgs), alpha is angle

465

c dens is electron density (cgs), alpha is angle

466

c between k and the magnetic field (radians)

466

c between k and the magnetic field (radians)

467

468

pi=4.0*atan(1.0)

468

pi=4.0*atan(1.0)

469

470

if(omega.eq.0.0) omega=1.0

470

if(omega.eq.0.0) omega=1.0

471

472

omegae=e*bfld*1.0e-4/em

472

omegae=e*bfld*1.0e-4/em

473

densmks=dens*1.0e6

473

densmks=dens*1.0e6

474

475

sum1=0.0

475

sum1=0.0

476

sum2=0.0

476

sum2=0.0

477

478

do i=1,nion

478

do i=1,nion

479

tr(i)=te/ti(i)

479

tr(i)=te/ti(i)

480

vti(i)=sqrt(bk*ti(i)/(1.67e-27*wi(i)))

480

vti(i)=sqrt(bk*ti(i)/(1.67e-27*wi(i)))

481

thetai(i)=(omega/ak)/(sqrt(2.0)*vti(i))

481

thetai(i)=(omega/ak)/(sqrt(2.0)*vti(i))

482

psii(i)=(vin(i)/ak)/(sqrt(2.0)*vti(i))

482

psii(i)=(vin(i)/ak)/(sqrt(2.0)*vti(i))

483

yi(i)=y_ion(thetai(i),psii(i))

483

yi(i)=y_ion(thetai(i),psii(i))

484

sum1=sum1+fi(i)*tr(i)*yi(i)

484

sum1=sum1+fi(i)*tr(i)*yi(i)

485

sum2=sum2+fi(i)*yi(i)

485

sum2=sum2+fi(i)*yi(i)

486

end do

486

end do

487

dl=ak**2*dlf*te/densmks

487

dl=ak**2*dlf*te/densmks

488

c write(*,fmt='("Before YE")')

489

c write(*,*) imode

490

c call exit

491

488

492

if(imode.eq.1.or.imode.eq.3) then

489

if(imode.eq.1.or.imode.eq.3) then

493

490

494

vte=sqrt(bk*te/em)

491

vte=sqrt(bk*te/em)

495

thetae=(omega/ak)/(sqrt(2.0)*vte)

492

thetae=(omega/ak)/(sqrt(2.0)*vte)

496

phi=(omegae/ak)/(sqrt(2.0)*vte)

493

phi=(omegae/ak)/(sqrt(2.0)*vte)

497

psie=(ven/ak)/(sqrt(2.0)*vte)

494

psie=(ven/ak)/(sqrt(2.0)*vte)

498

ye=y_electron(thetae,phi,psie,alpha)

495

ye=y_electron(thetae,phi,psie,alpha)

499

write(*,fmt='("AFTER YE")')

500

c call exit

501

496

502

else if(imode.eq.2) then

497

else if(imode.eq.2) then

503

c

498

c

504

c use Mike Sulzer's library here: alpha2 is angle off perp in degrees

499

c use Mike Sulzer's library here: alpha2 is angle off perp in degrees

505

c

500

c

506

501

507

freq=omega/(2.0*pi)

502

freq=omega/(2.0*pi)

508

alpha2=abs(pi/2.0-alpha)*180.0/pi

503

alpha2=abs(pi/2.0-alpha)*180.0/pi

509

c write(*,*) "ye: ", ye

504

c write(*,*) "ye: ", ye

510

call collision(densmks, te, freq, alpha2, ye)

505

call collision(densmks, te, freq, alpha2, ye)

511

c write(*,fmt='("AFTER COLLISION")')

506

c write(*,fmt='(" geobfield: time is before earliest model.")')

512

c call exit

513

ye=ye*omega+cmplx(0.0,1.0)

507

ye=ye*omega+cmplx(0.0,1.0)

514

508

515

end if

509

end if

516

yed=ye+cmplx(0.0,dl)

510

yed=ye+cmplx(0.0,dl)

517

511

518

p=(cabs(ye))**2*real(sum2)+cabs(sum1+cmplx(0.0,dl))**2*real(ye)

512

p=(cabs(ye))**2*real(sum2)+cabs(sum1+cmplx(0.0,dl))**2*real(ye)

519

p=p/(cabs(yed+sum1))**2

513

p=p/(cabs(yed+sum1))**2

520

spect1=p*2.0e0/(omega*pi)

514

spect1=p*2.0e0/(omega*pi)

521

write(*,*) "spect1:",spect1

515

522

return

516

return

523

end

517

end

524

518

525

subroutine acf2(wl, tau, te1, ti1, fi1, ven1, vin1, wi1,

519

subroutine acf2(wl, tau, te1, ti1, fi1, ven1, vin1, wi1, nion1,

526

& alpha1, dens1, bfld1, acf, ~~nion1~~)

520

& alpha1, dens1, bfld1, acf)

527

c

521

c

528

c computes autocorrelation function for given plasma parameters

522

c computes autocorrelation function for given plasma parameters

529

c by integrating real spectrum

523

c by integrating real spectrum

530

c tau in sec., alpha in radians, density in cgs, bfield in cgs

524

c tau in sec., alpha in radians, density in cgs, bfield in cgs

531

c scattering wavelength (wl) in meters

525

c scattering wavelength (wl) in meters

532

c

526

c

533

include 'fitter.h'

527

include 'fitter.h'

534

528

535

real wl,tau,te1,ti1(nion1),fi1(nion1),ven1,vin1(nion1),alpha1,

529

real wl,tau,te1,ti1(nion1),fi1(nion1),ven1,vin1(nion1),alpha1,

536

& dens1,bfld1,acf

530

& dens1,bfld1,acf

537

real pi

531

real pi

538

integer nion1

532

integer nion1

539

integer wi1(nion1)

533

integer wi1(nion1)

540

integer i,j,k,~~imode~~

534

integer i,j,k

541

common /mode/imode

542

c

535

c

543

c write(*,*) "INITIAL acf:",wl,tau,te1,ti1,fi1,ven1,vin1,wi1,alpha1

544

write(*,*) "INITIAL acf:",dens1, bfld1, acf, nion1

545

c write(*,fmt='("INIT")')

546

pi=4.0*atan(1.0)

536

pi=4.0*atan(1.0)

547

c

537

c

548

c copy arguments to common block

538

c copy arguments to common block

549

c

539

c

550

ak=2.0*pi/wl

540

ak=2.0*pi/wl

551

imode=2

541

imode=2

552

write(*,*) "imode:",imode

553

542

554

nion=nion1

543

nion=nion1

555

alpha=alpha1

544

alpha=alpha1

556

te=te1

545

te=te1

557

ven=ven1

546

ven=ven1

558

c write(*,fmt='("INIT2")')

559

do i=1,nion

547

do i=1,nion

560

c write(*,fmt='("INIT2.5")')

561

ti(i)=ti1(i)

548

ti(i)=ti1(i)

562

fi(i)=fi1(i)

549

fi(i)=fi1(i)

563

vin(i)=vin1(i)

550

vin(i)=vin1(i)

564

wi(i)=wi1(i)

551

wi(i)=wi1(i)

565

end do

552

end do

566

c write(*,fmt='("INIT3")')

567

dens=dens1

553

dens=dens1

568

bfld=bfld1

554

bfld=bfld1

569

555

570

c write(*,*) wl,alpha1,bfld1,dens

556

c write(*,*) wl,alpha1,bfld1,dens

571

c call exit

557

c call exit

572

c write(*,fmt='("Before Gauss")')

573

call gaussq(tau,acf)

574

558

575

write(*,*) "FINAL acf:",acf

559

call gaussq(tau,acf)

576

560

577

c write(*,fmt='("After Gauss")')

578

return

561

return

579

end

562

end

580

563

581

subroutine gaussq(tau,acf)

564

subroutine gaussq(tau,acf)

582

c

565

c

583

c Computes cosine or sine transform of given real function

566

c Computes cosine or sine transform of given real function

584

c Uses quadpack, tau in sec.

567

c Uses quadpack, tau in sec.

585

c

568

c

586

real a,abserr,epsabs,acf,tau,work,chebmo

569

real a,abserr,epsabs,acf,tau,work,chebmo

587

integer ier,integr,iwork,last,leniw,lenw,limit,limlst,

570

integer ier,integr,iwork,last,leniw,lenw,limit,limlst,

588

* lst,maxp1,neval

571

* lst,maxp1,neval

589

integer knorm,ksave,momcom,nrmom

572

integer knorm,ksave,momcom,nrmom

590

dimension iwork(300),work(1625),chebmo(61,25)

573

dimension iwork(300),work(1625),chebmo(61,25)

591

external spect1

574

external spect1

592

include 'fitter.h'

575

include 'fitter.h'

593

c

576

c

594

c write(*,*) "inside gaussq"

577

c write(*,*) "inside gaussq"

595

a = 0.0

578

a = 0.0

596

b = 2.5e4*ak ! upper integration limit open to debate (was 1.5, now 2.5)

579

b = 2.5e4*ak ! upper integration limit open to debate (was 1.5, now 2.5)

597

integr = 1

580

integr = 1

598

epsabs = 1.0e-4

581

epsabs = 1.0e-4

599

limlst = 50

582

limlst = 50

600

limit = 100

583

limit = 100

601

leniw = limit*2+limlst

584

leniw = limit*2+limlst

602

maxp1 = 61

585

maxp1 = 61

603

lenw = leniw*2+maxp1*25

586

lenw = leniw*2+maxp1*25

604

587

605

c write(*,*) leniw,lenw

588

c write(*,*) leniw,lenw

606

589

607

nrmom=0

590

nrmom=0

608

ksave=0

591

ksave=0

609

momcom=0

592

momcom=0

610

write(*,*) "Before qc25f:",acf,imode

593

611

c much faster, more robust

594

c much faster, more robust

612

c write(*,*) "acf_in: ",acf

595

c write(*,*) "acf_in: ",acf

613

call qc25f(spect1,a,b,tau,integr,nrmom,maxp1,ksave,acf,

596

call qc25f(spect1,a,b,tau,integr,nrmom,maxp1,ksave,acf,

614

& abserr,neval,resabs,resasc,momcom,chebmo)

597

& abserr,neval,resabs,resasc,momcom,chebmo)

615

c write(*,*) "acf_out: ",acf

598

c write(*,*) "acf_out: ",acf

616

c call qawf(spect1,a,tau,integr,epsabs,acf,abserr,neval,

599

c call qawf(spect1,a,tau,integr,epsabs,acf,abserr,neval,

617

c & ier,limlst,lst,leniw,maxp1,lenw,iwork,work)

600

c & ier,limlst,lst,leniw,maxp1,lenw,iwork,work)

618

write(*,*) "After qc25f:",acf

601

619

c call exit

602

620

return

603

return

621

end

604

end

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages

                   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)
 continue
 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)
 continue
 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
 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)
 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
-            c      write(*,fmt='("Before YE")')
-            c      write(*,*) imode
-            c      call exit
                   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)
-            	 write(*,fmt='("AFTER YE")')
-            c         call exit
                   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)
-            c         write(*,fmt='("AFTER COLLISION")')
+            c         write(*,fmt='(" geobfield: time is before earliest model.")')
-            c         call exit
                      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)
-                  write(*,*) "spect1:",spect1
                   return
                   end
-                  subroutine acf2(wl, tau, te1, ti1, fi1, ven1, vin1, wi1,
+                  subroutine acf2(wl, tau, te1, ti1, fi1, ven1, vin1, wi1, nion1,
-                 &     alpha1, dens1, bfld1, acf, nion1)
+                 &     alpha1, dens1, bfld1, acf)
             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)
-                  integer i,j,k,imode
+                  integer i,j,k
-                  common /mode/imode
             c
-            c      write(*,*) "INITIAL acf:",wl,tau,te1,ti1,fi1,ven1,vin1,wi1,alpha1
-                  write(*,*) "INITIAL acf:",dens1, bfld1, acf, nion1
-            c      write(*,fmt='("INIT")')
                   pi=4.0*atan(1.0)
             c
             c     copy arguments to common block
             c
                   ak=2.0*pi/wl
                   imode=2
-                  write(*,*) "imode:",imode
                   nion=nion1
                   alpha=alpha1
                   te=te1
                   ven=ven1
-            c      write(*,fmt='("INIT2")')
                   do i=1,nion
-            c      write(*,fmt='("INIT2.5")')
                     ti(i)=ti1(i)
                     fi(i)=fi1(i)
                     vin(i)=vin1(i)
                     wi(i)=wi1(i)
                   end do
-            c      write(*,fmt='("INIT3")')
                   dens=dens1
                   bfld=bfld1
             c      write(*,*) wl,alpha1,bfld1,dens
             c      call exit
-            c      write(*,fmt='("Before Gauss")')
-                  call gaussq(tau,acf)
-                  write(*,*) "FINAL acf:",acf
+                  call gaussq(tau,acf)
-            c      write(*,fmt='("After Gauss")')
                   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
-                  write(*,*) "Before qc25f:",acf,imode
             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)
-                  write(*,*) "After qc25f:",acf
-            c      call exit
                   return
                   end