          gam7=gam7*2.*(PH(I)/boltz/temp)
          gam7=gam7*(1.+heh*0.41336)
*     print *,gam7,' =FullWHM in rad/s'
          gamma=gamrad+gam7
        else
* Classical van der Waals broadening
          ALGH = 8.6735 + .7*ALOG10(THETA) + ALOG10(PH(I)) + .4*ALOGC6
          GAMH = EX10(ALGH)*1.E08
          GAM6 = GAMH*(1. + HEH*0.41336)*FDAMP
          GAMMA = GAMRAD + GAM6
        endif
*
        A(I) = GAMMA/(FYRAPI*DNUD(I))
*
    1 CONTINUE
*
      RETURN
      END
*
*
*
*
*
*
      FUNCTION gammln(xx) 
* returns ln gamma(xx) The Gamma function, Numerical recipes, 2nd ed, p 206
      REAL gammln,xx 
      INTEGER j 
      DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) 
      SAVE cof,stp 
      DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, 
     *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, 
     *-.5395239384953d-5,2.5066282746310005d0/ 
      x=xx 
      y=x 
      tmp=x+5.5d0 
      tmp=(x+0.5d0)*log(tmp)-tmp 
      ser=1.000000000190015d0 
      do 11 j=1,6 
        y=y+1.d0 
        ser=ser+cof(j)/y 
11    continue 
      gammln=tmp+log(stp*ser/x) 
      return 
      END 
C  (C) Copr. 1986-92 Numerical Recipes Software 1"n26@1kN. 
* see link ~/smallcalc/Numerical.recipes
*
*
*
*
*
      subroutine Anstee(iel,ion,xlb,elo,eion,sigma,velexp,levlo,levup,
     &                                                          recipe)
      save
* Get Hydrogen damping data for s-p and p-s transitions in neutral species
* ref Anstee S.D., O'Mara B.J 1995, MNRAS 276, 859-866
* table1 and table2 gives sigma(p,s) and velexp(p,s)
*
* INPUT parameters:
* iel   = atomic number
* ion   = ionization state (1 for neutral)
* xlb   = line wavelength in Angstroms
* elo   = excitation energy of the lower level (eV)
* eion  = ionization potential (eV)
* levlo = state designation of the lower energy level ( s, p, d etcetera)
* levup = state designation of the upper energy level ( s, p, d etcetera)
*
* OUTPUT parameters:
* sigma  = crossection in atomic units
* velexp = velocity exponent
* recipe = A for "Anstee" lines, U (Unsold) for all other lines
*
      implicit none
      real table1(13:30,10:30),table2(13:30,10:30),sigma,velexp
      real elo,eup,eion,xlb,nstarp,nstars
      integer i,j,ifirst,iel,ion
      character levlo,levup,string*132,check*1,recipe*1
*
      data ifirst / 1 /
*
      recipe='U'
      sigma=0.0
      velexp=0.0
      if(ion.ne.1) then
        return
* Anstee and O'Mara give data only for neutral species
      endif
*
*
      if(ifirst.eq.1) then
*
* Read s-p and p-s data from Anstee.dat for sigma and velexp
*
        open(612,file='data/Anstee.dat',status='old',form='formatted')
*
* First read data from table 1, not commentary lines
*
    1   continue
        read(612,1000) string
 1000   format(a132)
        read(string,1010) check
 1010   format(a1)
        if (check.ne.'*') then
          read(string,*) (table1(j,10),j=13,30)
          do i=11,30
            read(612,*) (table1(j,i),j=13,30)
          enddo
        else
          goto 1
        endif
*
* Read data from table 2, not commentary lines
*
    2   continue
        read(612,1000) string
        read(string,1010) check
        if (check.ne.'*') then
          read(string,*) (table2(j,10),j=13,30)
          do i=11,30
            read(612,*) (table2(j,i),j=13,30)
          enddo
        else
          goto 2
        endif
*
        close(612)
        ifirst=0
*
      endif
* 
*
      if ((levlo.eq.'s'.or.levlo.eq.'S').and.
     &    (levup.eq.'p'.or.levup.eq.'P')) then
        nstars=1./(sqrt((eion-elo)/13.60))
        eup=12398.426/xlb+elo
        nstarp=1./(sqrt((eion-eup)/13.60))
        if ((nstarp.ge.1.3.and.nstarp.le.3.0).and.
     &      (nstars.ge.1.0.and.nstars.le.3.0))then
          call interpolation(table1,table2,nstarp,nstars,
     &                                      sigma,velexp)
          recipe='A'
          return
        endif
      else if ((levlo.eq.'p'.or.levlo.eq.'P').and.
     &         (levup.eq.'s'.or.levup.eq.'S')) then
        nstarp=1./(sqrt((eion-elo)/13.60))
        eup=12398.426/xlb+elo
        nstars=1./(sqrt((eion-eup)/13.60))
        if ((nstarp.ge.1.3.and.nstarp.le.3.0).and.
     &      (nstars.ge.1.0.and.nstars.le.3.0)) then
          call interpolation(table1,table2,nstarp,nstars,        
     &                                     sigma,velexp)
          recipe='A'
          return
        endif
      endif
*
*
* Investigate whether the line was cpmputed by Paul Barklem
      call spec_lines(iel,ion,xlb,sigma,velexp,recipe)
      if(recipe.ne.'A') then
*
* OK, we have to stick to Unsolds recipe
*
        recipe='U'
      endif
*
      return
      end
*
************************************************************************
*
      subroutine interpolation(table1,table2,nstarp,nstars,sigma,velexp)
* interpolate in Anstee's tables to get cross section and velocity
*       exponent for s-p or p-s transitions
*
      implicit none
      real table1(13:30,10:30),table2(13:30,10:30),sigma,velexp
      real nstarp,nstars,fractp,fracts,a,b,c,d,x,y
      integer ip,is
*
      ip=int(10.*nstarp)
      is=int(10.*nstars)
      fractp=(10.*nstarp-float(ip))
      fracts=(10.*nstars-float(is))
* interpolate in cross section table, first in the p direction
      a=table1(ip,is)
      b=table1(ip+1,is)
      c=table1(ip,is+1)
      d=table1(ip+1,is+1)
      x=a*(1.-fractp)+b*fractp
      y=c*(1.-fractp)+d*fractp
* interpolate in the s direction
      sigma=x*(1.-fracts)+y*fracts
* interpolate in velocity exponent table, first in the p direction
      a=table2(ip,is)
      b=table2(ip+1,is)
      c=table2(ip,is+1)
      d=table2(ip+1,is+1)
      x=a*(1.-fractp)+b*fractp
      y=c*(1.-fractp)+d*fractp
* interpolate in the s direction
      velexp=x*(1.-fracts)+y*fracts
*
      return
      end
*
************************************************************************
*
      subroutine spec_lines(iel,ion,xlb,sigma,velexp,recipe)
*
* SPECIAL: a few p-d transitions computed by Paul Barklem 96.10.29 e-mail
*
      implicit none
      character*1 recipe
      integer iel,ion
      real xlb,sigma,velexp
*
      if(ion.ne.1) then
        return
* Peter Barklem gives data only for neutral species
      endif
*
*
      if(iel.eq.11) then
        if(abs(xlb-5682.65).le.0.05) then 
* Na 1 5682.65 AA test: nstarp= 2.116, nstard= 3.987
          recipe='A'
          velexp=    0.327
          sigma=    1969.
          return
        endif
        if(abs(xlb-5688.22).le.0.05) then 
* Na 1 5688.22 AA test: nstarp= 2.117, nstard= 3.987
          recipe='A'
          velexp=    0.327
          sigma=     1969.
          return
        endif
*
*
      else if(iel.eq.26) then
        if(abs(xlb-5679.03).le.0.05) then 
* Fe 1 5679.03 AA test: nstarp= 2.047, nstard= 3.576
          recipe='A'
          velexp=    0.306
          sigma=     1455.
          return
        endif
        if(abs(xlb-5686.54).le.0.05) then 
* Fe 1 5686.54 AA test: nstarp= 2.015, nstard= 3.410
          recipe='A'
          velexp=    0.303
          sigma=     1275.
          return
        endif
        if(abs(xlb-6170.52).le.0.05) then 
* Fe 1 6170.52 AA test: nstarp= 2.094, nstard= 3.526
          recipe='A'
          velexp=    0.309
          sigma=     1406.
          return
        endif
        if(abs(xlb-7751.12).le.0.05) then 
* Fe 1 7751.12 AA test: nstarp= 2.163, nstard= 3.225
          recipe='A'
          velexp=    0.294
          sigma=     1042.
          return
        endif
        if(abs(xlb-7807.92).le.0.05) then 
* Fe 1 7807.92 AA test: nstarp= 2.163, nstard= 3.211
          recipe='A'
          velexp=    0.284
          sigma=     1022.
          return
        endif
*
*
      else if(iel.eq.28) then
        if(abs(xlb-5682.21).le.0.05) then 
* Ni 1 5682.21 AA test: nstarp= 1.963, nstard= 3.176
          recipe='A'
          velexp=    0.280
          sigma=     1028.
          return
        endif
        if(abs(xlb-5694.99).le.0.05) then 
* Ni 1 5694.99 AA test: nstarp= 1.958, nstard= 3.152
          recipe='A'
          velexp=    0.263
          sigma=     1000.
          return
        endif
        if(abs(xlb-6130.14).le.0.05) then 
* Ni 1 6130.14 AA test: nstarp= 2.009, nstard= 3.118
          recipe='A'
          velexp=    0.277
          sigma=      960.
          return
        endif
        if(abs(xlb-6175.37).le.0.05) then 
* Ni 1 6175.37 AA test: nstarp= 1.958, nstard= 2.973
          recipe='A'
          velexp=    0.269
          sigma=      792.
          return
        endif
        if(abs(xlb-6176.82).le.0.05) then 
* Ni 1 6176.82 AA test: nstarp= 1.958, nstard= 2.972
          recipe='A'
          velexp=    0.268
          sigma=      791.
          return
        endif
      endif
*
      return
      end
*
************************************************************************





*
*	This is data from S.D. Anstee and B.J. O'Mara: Mon. Not. R. Aston. Soc. 276, 859-866 (1995)
*	"Width cross-sections for collisional broadening of s-p and p-s transitions by atomic hydrogen"
*       Multiply Table 1 data by the Bohr radius area=a0**2=2.8002e-17 to get units of cm-2
*
*	Table 1: Line broadening cross-section for relative collision velocity 10^4 m/s.
*
126. 140. 165. 202. 247. 299. 346. 383. 435. 491. 553. 617. 685. 769. 838. 925. 1011. 1082.
140. 150. 162. 183. 218. 273. 327. 385. 440. 501. 557. 620. 701. 764. 838. 923. 1025. 1085.
154. 167. 175. 192. 216. 251. 299. 357. 423. 487. 549. 617. 684. 759. 834. 910. 1014. 1064.
166. 180. 192. 206. 226. 253. 291. 339. 397. 459. 532. 600. 676. 755. 832. 896. 1002. 1055.
208. 194. 207. 223. 242. 265. 296. 335. 384. 445. 511. 583. 656. 726. 817. 889. 988. 1044.
262. 254. 220. 239. 261. 283. 310. 344. 388. 442. 496. 568. 635. 725. 791. 890. 970. 1036.
311. 306. 299. 251. 280. 304. 330. 361. 396. 443. 500. 563. 630. 704. 796. 880. 951. 1033.
358. 359. 350. 338. 293. 323. 352. 381. 416. 455. 511. 566. 635. 706. 780. 859. 946. 1039.
411. 409. 405. 392. 370. 340. 375. 406. 439. 478. 525. 580. 644. 714. 790. 873. 961. 1050.
462. 463. 459. 450. 443. 400. 394. 432. 467. 501. 546. 595. 650. 711. 786. 873. 963. 1050.
522. 525. 529. 524. 516. 518. 438. 454. 495. 532. 565. 621. 671. 741. 813. 874. 951. 1034.
589. 593. 590. 583. 579. 568. 565. 483. 517. 560. 600. 644. 691. 752. 821. 904. 978. 1048.
658. 655. 666. 657. 649. 653. 649. 587. 549. 592. 674. 674. 728. 782. 833. 902. 992. 1084.
738. 742. 747. 725. 721. 729. 699. 730. 626. 622. 668. 721. 765. 809. 887. 938. 1001. 1109.
838. 838. 810. 809. 790. 800. 769. 815. 757. 679. 704. 755. 806. 854. 901. 974. 1034. 1105.
942. 946. 925. 901. 918. 895. 919. 897. 933. 890. 785. 797. 859. 908. 976. 1020. 1115. 1173.
1059. 1061. 1056. 1061. 1074. 1031. 1036. 1036. 993. 1038. 932. 852. 878. 943. 1003. 1074. 1131. 1200.
1069. 1076. 1083. 1095. 1102. 1091. 1126. 1156. 1103. 1149. 1157. 1036. 972. 1007. 1064. 1124. 1209. 1283.
1338. 1350. 1356. 1354. 1324. 1301. 1312. 1318. 1257. 1239. 1297. 1233. 1089. 1059. 1106. 1180. 1218. 1317.
1409. 1398. 1367. 1336. 1313. 1313. 1409. 1354. 1317. 1287. 1353. 1386. 1279. 1158. 1141. 1188. 1260. 1335.
1328. 1332. 1342. 1369. 1405. 1451. 1502. 1524. 1506. 1477. 1522. 1594. 1572. 1436. 1328. 1325. 1382. 1446.
*
*
*
*	Table 2: Velocity exponents alfa corresponding to the cross-section of table 1.
*
.268 .269 .335 .377 .327 .286 .273 .270 .271 .268 .267 .264 .264 .264 .261 .256 .248 .245
.261 .256 .254 .282 .327 .355 .321 .293 .287 .271 .267 .273 .270 .270 .268 .268 .264 .263
.266 .264 .257 .252 .267 .289 .325 .339 .319 .301 .292 .284 .281 .281 .277 .282 .276 .274
.262 .274 .258 .251 .247 .254 .273 .291 .316 .322 .320 .302 .294 .290 .287 .292 .283 .277
.322 .275 .264 .259 .250 .245 .273 .255 .271 .284 .294 .308 .296 .299 .288 .289 .282 .278
.267 .300 .260 .268 .254 .242 .243 .242 .239 .246 .267 .277 .280 .290 .282 .281 .274 .271
.259 .274 .275 .252 .265 .248 .249 .237 .283 .236 .247 .254 .254 .271 .268 .267 .258 .262
.260 .255 .268 .268 .268 .264 .248 .239 .229 .240 .236 .234 .238 .244 .252 .251 .244 .255
.255 .255 .244 .247 .317 .246 .255 .244 .237 .231 .227 .231 .235 .232 .235 .241 .237 .245
.256 .254 .254 .249 .227 .319 .253 .253 .240 .237 .238 .233 .231 .230 .228 .234 .227 .241
.257 .254 .252 .235 .253 .240 .284 .251 .246 .241 .235 .228 .222 .225 .225 .219 .228 .233
.244 .240 .245 .238 .248 .230 .283 .252 .244 .244 .238 .235 .234 .236 .228 .224 .225 .231
.224 .241 .244 .237 .237 .249 .219 .324 .239 .245 .242 .242 .232 .233 .221 .227 .231 .218
.241 .245 .249 .239 .243 .250 .217 .254 .308 .237 .247 .244 .234 .228 .233 .224 .227 .226
.243 .243 .232 .227 .235 .253 .227 .220 .320 .270 .243 .252 .248 .238 .234 .241 .225 .227
.225 .226 .234 .230 .226 .233 .249 .225 .216 .300 .286 .237 .240 .247 .243 .234 .231 .238
.268 .260 .247 .238 .233 .241 .254 .248 .207 .227 .315 .260 .226 .237 .240 .239 .239 .240
.248 .246 .238 .226 .213 .221 .226 .226 .204 .194 .248 .316 .234 .216 .236 .233 .221 .230
.200 .202 .198 .194 .206 .207 .227 .224 .207 .185 .198 .275 .315 .233 .229 .231 .233 .236
.202 .209 .221 .226 .230 .245 .202 .257 .246 .225 .215 .246 .320 .321 .244 .239 .251 .253
.246 .248 .255 .265 .274 .285 .292 .284 .273 .250 .225 .239 .295 .352 .320 .258 .260 .269

