| 1 | subroutine geteii(z_in, ion_in, temp_in, rate)
|
|---|
| 2 |
|
|---|
| 3 | c================
|
|---|
| 4 | c
|
|---|
| 5 | c For a given ionisation stage (ion_in) of an element (z_in)
|
|---|
| 6 | c this code returns the electron impact ionisation (EII) rate
|
|---|
| 7 | c coefficient (rate) at temperature temp_in. EII data comes from
|
|---|
| 8 | c Dere, 2007, A&A, 466, 771 and this codes calls a cubic spline
|
|---|
| 9 | c routine to return rate at the given temp_in.
|
|---|
| 10 | c
|
|---|
| 11 | c Inputs:
|
|---|
| 12 | c z_in = elements (1-30, i.e., H-Zn)
|
|---|
| 13 | c ion_in = charge state of ion AFTER ionisation (1-30)
|
|---|
| 14 | c temp_in = temperature at which rate coefficient is desired (K)
|
|---|
| 15 | c
|
|---|
| 16 | c Outputs:
|
|---|
| 17 | c rate = EII rate coefficient at temp_in (cm^-3 s^-1)
|
|---|
| 18 | c
|
|---|
| 19 | c
|
|---|
| 20 | c Paul Bryans
|
|---|
| 21 | c Columbia University
|
|---|
| 22 | c 16 September 2008
|
|---|
| 23 | c
|
|---|
| 24 | c================
|
|---|
| 25 |
|
|---|
| 26 | implicit none
|
|---|
| 27 |
|
|---|
| 28 | real*8 temp(300), eii(300), temp_in, rate
|
|---|
| 29 | real*8 logtemp(300), logeii(300), logtemp_in, logeii_in
|
|---|
| 30 | real*8 d2logeii(300)
|
|---|
| 31 | logical lspl
|
|---|
| 32 | character eiifile*120, element(30)*2
|
|---|
| 33 | integer i, z_in, ion_in, npts
|
|---|
| 34 |
|
|---|
| 35 | data(element(i),i=1,30)/'h','he','li','be','b','c','n','o','f',
|
|---|
| 36 | & 'ne','na','mg','al','si','p','s','cl',
|
|---|
| 37 | & 'ar','k','ca','sc','ti','v','cr','mn',
|
|---|
| 38 | & 'fe','co','ni','cu','zn'/
|
|---|
| 39 |
|
|---|
| 40 |
|
|---|
| 41 |
|
|---|
| 42 | logtemp_in=dlog10(temp_in)
|
|---|
| 43 | lspl=.true.
|
|---|
| 44 |
|
|---|
| 45 |
|
|---|
| 46 | if (ion_in.gt.9) then
|
|---|
| 47 | eiifile='/path to files here/'//
|
|---|
| 48 | & element(z_in)(1:len_trim(element(z_in)))//'_'//
|
|---|
| 49 | & char(mod(ion_in/10,10)+48)//char(mod(ion_in,10)+48)
|
|---|
| 50 | & //'_chianti.ionizdat'
|
|---|
| 51 | else
|
|---|
| 52 | eiifile='/path to files here/'//
|
|---|
| 53 | & element(z_in)(1:len_trim(element(z_in)))//'_'//
|
|---|
| 54 | & char(mod(ion_in,10)+48)//'_chianti.ionizdat'
|
|---|
| 55 | endif
|
|---|
| 56 |
|
|---|
| 57 | write(6,*)eiifile
|
|---|
| 58 |
|
|---|
| 59 | call read_ionrecdat(eiifile,temp,eii,npts)
|
|---|
| 60 |
|
|---|
| 61 |
|
|---|
| 62 | if(lspl)then
|
|---|
| 63 |
|
|---|
| 64 |
|
|---|
| 65 | do i=1,npts
|
|---|
| 66 | logtemp(i)=dlog10(temp(i))
|
|---|
| 67 | logeii(i)=dlog10(eii(i))
|
|---|
| 68 | enddo
|
|---|
| 69 |
|
|---|
| 70 | call spline(npts, logtemp, logeii, d2logeii)
|
|---|
| 71 | call spline_int(npts, logtemp, logeii, d2logeii,
|
|---|
| 72 | & logtemp_in, logeii_in)
|
|---|
| 73 |
|
|---|
| 74 | else
|
|---|
| 75 |
|
|---|
| 76 | c linear interpolation not implemented
|
|---|
| 77 |
|
|---|
| 78 | endif
|
|---|
| 79 |
|
|---|
| 80 | if(logeii_in.eq.0)then
|
|---|
| 81 | rate=0d0
|
|---|
| 82 | else
|
|---|
| 83 | rate =10d0**(logeii_in)
|
|---|
| 84 | endif
|
|---|
| 85 |
|
|---|
| 86 |
|
|---|
| 87 | return
|
|---|
| 88 |
|
|---|
| 89 |
|
|---|
| 90 | end
|
|---|
| 91 |
|
|---|
| 92 | c================
|
|---|
| 93 | c================
|
|---|
| 94 |
|
|---|
| 95 |
|
|---|
| 96 | subroutine read_ionrecdat(file,temp,rate,npts)
|
|---|
| 97 |
|
|---|
| 98 | c================
|
|---|
| 99 | c
|
|---|
| 100 | c Returns temperature and rate coefficient values from EII
|
|---|
| 101 | c tables of Dere, 2007, A&A, 466, 771.
|
|---|
| 102 | c
|
|---|
| 103 | c Paul Bryans
|
|---|
| 104 | c Columbia University
|
|---|
| 105 | c 16 September 2008
|
|---|
| 106 | c
|
|---|
| 107 | c================
|
|---|
| 108 |
|
|---|
| 109 | implicit none
|
|---|
| 110 |
|
|---|
| 111 | real*8 temp(300), rate(300)
|
|---|
| 112 | character file*120
|
|---|
| 113 | integer i, npts, iunt12
|
|---|
| 114 | logical open12
|
|---|
| 115 | parameter (iunt12=12)
|
|---|
| 116 |
|
|---|
| 117 |
|
|---|
| 118 | if (open12) then
|
|---|
| 119 | close(12)
|
|---|
| 120 | open12=.false.
|
|---|
| 121 | endif
|
|---|
| 122 | open( unit=iunt12 , file=file , status='unknown')
|
|---|
| 123 | open12=.true.
|
|---|
| 124 |
|
|---|
| 125 | i=0
|
|---|
| 126 | 10 continue
|
|---|
| 127 | i=i+1
|
|---|
| 128 | read(iunt12,1000) temp(i), rate(i)
|
|---|
| 129 | if(temp(i).gt.0) goto 10
|
|---|
| 130 |
|
|---|
| 131 | npts=i-1
|
|---|
| 132 |
|
|---|
| 133 | return
|
|---|
| 134 |
|
|---|
| 135 | 1000 format(2e12.3)
|
|---|
| 136 |
|
|---|
| 137 | end
|
|---|
| 138 |
|
|---|
| 139 | c================
|
|---|
| 140 | c================
|
|---|
| 141 |
|
|---|
| 142 |
|
|---|
| 143 | subroutine spline(npts, xin, yin, d2y)
|
|---|
| 144 |
|
|---|
| 145 | c================
|
|---|
| 146 | c
|
|---|
| 147 | c For an array of xin and yin values, this routine
|
|---|
| 148 | c calculates the second derivatives of y (d2y) for
|
|---|
| 149 | c cubic spline interpolation with d2y set to zero
|
|---|
| 150 | c at both end points.
|
|---|
| 151 | c
|
|---|
| 152 | c Paul Bryans
|
|---|
| 153 | c Columbia University
|
|---|
| 154 | c 16 September 2008
|
|---|
| 155 | c
|
|---|
| 156 | c================
|
|---|
| 157 |
|
|---|
| 158 |
|
|---|
| 159 | implicit none
|
|---|
| 160 |
|
|---|
| 161 | integer npts, i, nmax
|
|---|
| 162 | parameter (nmax=300)
|
|---|
| 163 | real*8 xin(nmax), yin(nmax), d2y(nmax)
|
|---|
| 164 | real*8 a(nmax), b(nmax), c(nmax), d(nmax)
|
|---|
| 165 | real*8 bb, tmp(nmax)
|
|---|
| 166 |
|
|---|
| 167 |
|
|---|
| 168 | c================
|
|---|
| 169 | c
|
|---|
| 170 | c natural cubic spline with second derivatives at end points set to zero
|
|---|
| 171 | c
|
|---|
| 172 | c================
|
|---|
| 173 |
|
|---|
| 174 | d2y(1) = 0d0
|
|---|
| 175 | d2y(npts) = 0d0
|
|---|
| 176 |
|
|---|
| 177 |
|
|---|
| 178 | c================
|
|---|
| 179 | c
|
|---|
| 180 | c put tridiagonal system of second derivatives in the form,
|
|---|
| 181 | c a(i)*d2y(i-1) + b(i)*d2y(i) + c(i)*d2y(i+1) = d(i)
|
|---|
| 182 | c
|
|---|
| 183 | c================
|
|---|
| 184 |
|
|---|
| 185 |
|
|---|
| 186 | do i=2,npts-1
|
|---|
| 187 | a(i) = (xin(i)-xin(i-1))/6d0
|
|---|
| 188 | b(i) = (xin(i+1)-xin(i-1))/3d0
|
|---|
| 189 | c(i) = (xin(i+1)-xin(i))/6d0
|
|---|
| 190 | d(i) = (yin(i+1)-yin(i))/(xin(i+1)-xin(i)) -
|
|---|
| 191 | & (yin(i)-yin(i-1))/(xin(i)-xin(i-1))
|
|---|
| 192 | enddo
|
|---|
| 193 |
|
|---|
| 194 |
|
|---|
| 195 |
|
|---|
| 196 | c================
|
|---|
| 197 | c
|
|---|
| 198 | c forward substitution of the tridiagonal system
|
|---|
| 199 | c
|
|---|
| 200 | c================
|
|---|
| 201 |
|
|---|
| 202 | bb = b(2)
|
|---|
| 203 | d2y(2) = d(2)/bb
|
|---|
| 204 | do i=3,npts-1
|
|---|
| 205 | tmp(i) = c(i-1)/bb
|
|---|
| 206 | bb = b(i) - a(i)*tmp(i)
|
|---|
| 207 | d2y(i) = (d(i)-a(i)*d2y(i-1))/bb
|
|---|
| 208 | enddo
|
|---|
| 209 |
|
|---|
| 210 |
|
|---|
| 211 | c================
|
|---|
| 212 | c
|
|---|
| 213 | c backward substitution of the tridiagonal system
|
|---|
| 214 | c
|
|---|
| 215 | c================
|
|---|
| 216 |
|
|---|
| 217 | do i=npts-2,2,-1
|
|---|
| 218 | d2y(i)=d2y(i)-c(i)*d2y(i+1)
|
|---|
| 219 | enddo
|
|---|
| 220 |
|
|---|
| 221 | return
|
|---|
| 222 |
|
|---|
| 223 | end
|
|---|
| 224 |
|
|---|
| 225 |
|
|---|
| 226 |
|
|---|
| 227 |
|
|---|
| 228 | c================
|
|---|
| 229 | c================
|
|---|
| 230 |
|
|---|
| 231 | subroutine spline_int(npts, xin, yin, d2y, xout, yout)
|
|---|
| 232 |
|
|---|
| 233 | c================
|
|---|
| 234 | c
|
|---|
| 235 | c For an array of xin and yin values, d2y is the
|
|---|
| 236 | c second derivative of yin calculated by spline
|
|---|
| 237 | c routine. This routine returns yout for a given
|
|---|
| 238 | c xout by interpolating.
|
|---|
| 239 | c
|
|---|
| 240 | c Paul Bryans
|
|---|
| 241 | c Columbia University
|
|---|
| 242 | c 16 September 2008
|
|---|
| 243 | c
|
|---|
| 244 | c================
|
|---|
| 245 |
|
|---|
| 246 |
|
|---|
| 247 | implicit none
|
|---|
| 248 |
|
|---|
| 249 | integer npts, i, index, nmax
|
|---|
| 250 | parameter (nmax=300)
|
|---|
| 251 | real*8 xin(nmax), yin(nmax), xout, yout
|
|---|
| 252 | real*8 d2y(nmax), a, b, c, d
|
|---|
| 253 |
|
|---|
| 254 |
|
|---|
| 255 | c================
|
|---|
| 256 | c
|
|---|
| 257 | c extrapolation
|
|---|
| 258 | c
|
|---|
| 259 | c================
|
|---|
| 260 |
|
|---|
| 261 | if (xout.lt.xin(1)) then
|
|---|
| 262 | write(6,*) 'warning: extrapolation '
|
|---|
| 263 | index=1
|
|---|
| 264 | a=1d0
|
|---|
| 265 | b=0d0
|
|---|
| 266 | c=xout-xin(1)
|
|---|
| 267 | d=0d0
|
|---|
| 268 | elseif (xout.gt.xin(npts)) then
|
|---|
| 269 | write(6,*) 'warning: extrapolation '
|
|---|
| 270 | index=npts-1
|
|---|
| 271 | a=0d0
|
|---|
| 272 | b=1d0
|
|---|
| 273 | c=0d0
|
|---|
| 274 | d=xout-xin(npts)
|
|---|
| 275 |
|
|---|
| 276 |
|
|---|
| 277 | c================
|
|---|
| 278 | c
|
|---|
| 279 | c find position of xout in xin array for interpolation
|
|---|
| 280 | c
|
|---|
| 281 | c================
|
|---|
| 282 |
|
|---|
| 283 | else
|
|---|
| 284 | do i=npts,2,-1
|
|---|
| 285 | if(xout.le.xin(i)) index=i-1
|
|---|
| 286 | enddo
|
|---|
| 287 |
|
|---|
| 288 |
|
|---|
| 289 | c================
|
|---|
| 290 | c
|
|---|
| 291 | c set up interpolation parameters,
|
|---|
| 292 | c yout = a*yin(i) + b*yin(i+1) + c*d2y(i) + c*d2y(i+1)
|
|---|
| 293 | c
|
|---|
| 294 | c================
|
|---|
| 295 |
|
|---|
| 296 | a = (xin(index+1)-xout)/(xin(index+1)-xin(index))
|
|---|
| 297 | b = 1d0-a
|
|---|
| 298 | c = ((a**3d0-a)*(xin(index+1)-xin(index))**2d0)/6d0
|
|---|
| 299 | d = ((b**3d0-b)*(xin(index+1)-xin(index))**2d0)/6d0
|
|---|
| 300 |
|
|---|
| 301 | endif
|
|---|
| 302 |
|
|---|
| 303 | yout = a*yin(index) + b*yin(index+1) + c*d2y(index)
|
|---|
| 304 | & + d*d2y(index+1)
|
|---|
| 305 |
|
|---|
| 306 |
|
|---|
| 307 | return
|
|---|
| 308 |
|
|---|
| 309 | end
|
|---|
| 310 |
|
|---|
