SUBROUTINE COLHE(T,COL) C ======================= C C Helium (both neutral and ionized) collision rates C C Meaning of ICOL: for all kinds of transitions and for both HeI C and HeII: C ICOL = 0 - approximate expressions taken from Mihalas, Heasley C and Auer, NCAR-TN/STR-104 (1975) - for He I and II C C New expression for He II collisional ionization from Klaus Werner C (ICOL >= 1) (also originally from Mihalas) C C For He I bound-bound transitions, the following standard C possibilities are also available: C C ICOL = 1, 2, or 3 - much more accurate Storey's rates, C subroutine written by D.G.Hummer (COLLHE). C This procedure can be used only for transitions C between states with n = 1, 2, 3, 4. C ICOL = 1 - means that a given transition is a transition C between non-averaged (l,s) states. In this case, C labeling of the He I energy levels must agree C with that given in subroutine COLLHE, ie. states C have to be labeled sequentially in order of C increasing frequency. C ICOL = 2 - means that a given transition is a transition between C a non-averaged (l,s) lower state and averaged upper C state. C ICOL = 3 - means that a given transition is a transition between C two averaged states. C Note: C The program allows only two standard possibilities of C constructing averaged levels of He I: C i) all states within given principal quantum number n (>1) are C lumped together C ii) all siglet states for given n, and all triplet states for C given n are lumped together separately (there are thus two C explicit levels for a given n). C If the user wants to use another averaging, he had to take care C of appropriate averaged collisional rates himself (by updating C subroutine CSPEC) C C C ICOL < 0 - non-standard, user supplied formula C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'ATOMIC.FOR' DIMENSION COL(MTRANS) DIMENSION FHE1(16),G0(3),G1(3),G2(3),G3(3),A(6,10) PARAMETER (EXPIA1=-0.57721566,EXPIA2=0.99999193, * EXPIA3=-0.24991055,EXPIA4=0.05519968, * EXPIA5=-0.00976004,EXPIA6=0.00107857, * EXPIB1=0.2677734343,EXPIB2=8.6347608925, * EXPIB3=18.059016973,EXPIB4=8.5733287401, * EXPIC1=3.9584969228,EXPIC2=21.0996530827, * EXPIC3=25.6329561486,EXPIC4=9.5733223454) DATA FHE1/0.,2.75D-1,7.29D-2,2.96D-2,1.48 D-2,8.5D-3,5.3D-3, * 3.5D-3,2.5D-3,1.8D-3,1.5D-3,1.2D-3,9.4D-4,7.5D-4, * 6.1D-4,5.3D-4/ DATA G0/ 7.3399521D-2, 1.7252867, 8.6335087 /, * G1/-1.4592763D-7, 2.0944117D-6, 2.7575544D-5 /, * G2/ 7.6621299D5, 5.4254879D6, 6.6395519D6 /, * G3/ 2.3775439D2, 2.2177891D3, 5.20725D3 / DATA ((A(I,J),J=1,10),I=1,6) / * -8.5931587 , 85.014091 , 923.64099, 2018.6470, 1551.5061 , * -2327.4819 ,-10701.481 ,-27619.789,-41099.602,-61599.023 , * 9.3868790 ,-78.834488 ,-969.18451,-2243.1768,-2059.9768 , * 1546.7107 , 9834.3447 , 27067.436, 41421.254, 63594.133 , * -4.0027571 , 28.360615 , 401.23965, 983.83374, 1051.4103 , * -204.82320 ,-3335.4211 ,-10100.119,-15863.257,-24949.125 , * 0.83941799 ,-4.7963457 ,-81.122566,-209.86169,-251.30855 , * -43.175175 , 530.37292 , 1826.1049, 2941.6460, 4740.8364 , * -.86396709E-01,0.37385577 , 8.0078983, 21.757591, 28.375637 , * 11.890312 ,-39.536087 ,-161.52513,-266.86011,-440.88257 , * 0.34853835E-02,-.10401310E-01,-.30957383,-.87988985,-1.2254572 , * -.72724497 , 1.0879648 , 5.6239786, 9.5323009, 16.150818 / SAVE FHE1,G0,G1,G2,G3 C HKT=HK/T TK=HKT/H SRT=SQRT(T) t0=t CT=5.465D-11*SRT CT1=5.4499487/T/SRT C C -------------- C Neutral helium C -------------- C IF(IELHE1.EQ.0) GO TO 60 ICALL=0 N0I=NFIRST(IELHE1) N1I=NLAST(IELHE1) NKI=NNEXT(IELHE1) N0Q=NQUANT(NLAST(IELHE1))+1 N1Q=ICUP(IELHE1) DO 50 II=N0I,N1I IT=ITRA(II,NKI) IF(IT.EQ.0) GO TO 10 C C ******** Collisional ionization C IC=ICOL(IT) C1=OSC0(IT) C2=CPAR(IT) U0=ENION(II)*TK IF(IC.GE.0) THEN U1=U0+0.27 U2=(U0+3.43)/(U0+1.43)**3 IF(U0.LE.UN) THEN EXPIU0=-LOG(U0)+EXPIA1+U0*(EXPIA2+U0*(EXPIA3+U0*(EXPIA4+ * U0*(EXPIA5+U0*EXPIA6)))) ELSE EXPIU0=EXP(-U0)*((EXPIB1+U0*(EXPIB2+U0*(EXPIB3+ * U0*(EXPIB4+U0))))/(EXPIC1+U0*(EXPIC2+ * U0*(EXPIC3+U0*(EXPIC4+U0)))))/U0 END IF IF(U1.LE.UN) THEN EXPIU1=-LOG(U1)+EXPIA1+U1*(EXPIA2+U1*(EXPIA3+U1*(EXPIA4+ * U1*(EXPIA5+U1*EXPIA6)))) ELSE EXPIU1=EXP(-U1)*((EXPIB1+U1*(EXPIB2+U1*(EXPIB3+ * U1*(EXPIB4+U1))))/(EXPIC1+U1*(EXPIC2+ * U1*(EXPIC3+U1*(EXPIC4+U1)))))/U1 END IF COL(IT)=CT*C1*U0*(EXPIU0-U0*(0.728*EXPIU1/U1+ * 0.189*EXP(-U0)*U2)) ELSE CALL CSPEC(II,NKI,IC,C1,C2,U0,T,COL(IT)) END IF 10 IF(II.GE.N1I) GO TO 30 C C ********* Collisional excitation C DO 20 JJ=II+1,N1I ICT=ITRA(II,JJ) IF(ICT.EQ.0) GO TO 20 IC=ICOL(ICT) C1=OSC0(ICT) C2=CPAR(ICT) U0=FR0(ICT)*HKT IF(IC.EQ.0) THEN C C *** ICOL = 0 Formula used by Mihalas, Heasley, and Auer C IF(U0.LE.UN) THEN EX=-LOG(U0)+EXPIA1+U0*(EXPIA2+U0*(EXPIA3+U0*(EXPIA4+ * U0*(EXPIA5+U0*EXPIA6)))) ELSE EX=EXP(-U0)*((EXPIB1+U0*(EXPIB2+U0*(EXPIB3+ * U0*(EXPIB4+U0))))/(EXPIC1+U0*(EXPIC2+ * U0*(EXPIC3+U0*(EXPIC4+U0)))))/U0 END IF IF(II.EQ.N0I) THEN C C excitation from the ground state C COL(ICT)=CT1*EX/U0*C1 ELSE C C transitions between excited states C U1=U0+0.2 IF(U1.LE.UN) THEN EXPIU1=-LOG(U1)+EXPIA1+U1*(EXPIA2+U1*(EXPIA3+U1*(EXPIA4+ * U1*(EXPIA5+U1*EXPIA6)))) ELSE EXPIU1=EXP(-U1)*((EXPIB1+U1*(EXPIB2+U1*(EXPIB3+ * U1*(EXPIB4+U1))))/(EXPIC1+U1*(EXPIC2+ * U1*(EXPIC3+U1*(EXPIC4+U1)))))/U1 END IF COL(ICT)=CT1/U0*(EX-U0/U1*0.81873*EXPIU1)*C1 END IF ELSE IF(IC.EQ.1) THEN C C *** ICOL = 1 Storey - Hummer collisional rates between C non-averaged states C (Note: procedure COLLHE, which calculates all rates, C is called only once) C IF(ICALL.EQ.0) CALL COLLHE(T,COLHE1) ICALL=1 COL(ICT)=COLHE1(II-N0I+1,JJ-N0I+1) ELSE IF(IC.EQ.2.OR.IC.EQ.3) THEN C C *** ICOL = 2 or 3 Storey - Hummer collisional rates between C averaged states C IF(ICALL.EQ.0) CALL COLLHE(T,COLHE1) ICALL=1 COL(ICT)=CHEAV(II,JJ,IC) ELSE IF(IC.LT.0) THEN C C Non-standard, user supplied formula C CALL CSPEC(II,JJ,IC,C1,CPAR(ICT),U0,T,COL(ICT)) END IF 20 CONTINUE C C collisional excitations from level II to higher, non-explicit C levels are lumped into the collisional ionization rate C (the so-called modified collision ionization rate); C the individual rates are calculated by expressions used by C Mihalas, Heasley, and Auer C 30 IF(N1Q.EQ.0.OR.IT.EQ.0) GO TO 50 I=NQUANT(II) REL=G(II)/2./I/I DO 40 J=N0Q,N1Q XJ=J U0=(ENION(II)-EH/XJ/XJ)*TK IF(I.EQ.1) THEN GAM=0. C1=FHE1(J) ELSE C1=OSH(I,J)*REL U1=U0+0.2 IF(U1.LE.UN) THEN EXPIU1=-LOG(U1)+EXPIA1+U1*(EXPIA2+U1*(EXPIA3+U1*(EXPIA4+ * U1*(EXPIA5+U1*EXPIA6)))) ELSE EXPIU1=EXP(-U1)*((EXPIB1+U1*(EXPIB2+U1*(EXPIB3+ * U1*(EXPIB4+U1))))/(EXPIC1+U1*(EXPIC2+ * U1*(EXPIC3+U1*(EXPIC4+U1)))))/U1 END IF GAM=U0/U1*0.81873*EXPIU1 END IF IF(U0.LE.UN) THEN EXPIU0=-LOG(U0)+EXPIA1+U0*(EXPIA2+U0*(EXPIA3+U0*(EXPIA4+ * U0*(EXPIA5+U0*EXPIA6)))) ELSE EXPIU0=EXP(-U0)*((EXPIB1+U0*(EXPIB2+U0*(EXPIB3+ * U0*(EXPIB4+U0))))/(EXPIC1+U0*(EXPIC2+ * U0*(EXPIC3+U0*(EXPIC4+U0)))))/U0 END IF COL(IT)=COL(IT)+CT1/U0*C1*(EXPIU0-GAM) 40 CONTINUE 50 CONTINUE C C -------------- C Ionized helium C -------------- C 60 IF(IELHE2.EQ.0) RETURN N0I=NFIRST(IELHE2) N1I=NLAST(IELHE2) NKI=NNEXT(IELHE2) N0Q=NQUANT(NLAST(IELHE2))+1 N1Q=ICUP(IELHE2) X=LOG10(T) X2=X*X X3=X2*X X4=X3*X X5=X4*X CT2=3.7036489/T/SRT C DO 200 II=N0I,N1I I=II-N0I+1 IT=ITRA(II,NKI) IF(IT.EQ.0) GO TO 100 C C ********* Collisional ionization C c for high temperature, use XSTAR formulae C if(t0.gt.1.e5) then rno=16. izc=2 call irc(i,t0,izc,rno,cs) col(it)=cs go to 100 end if C IC=ICOL(IT) U0=FR0(IT)*HKT IF(IC.EQ.0) THEN IF(I.LE.3) THEN GAM=G0(I)-G1(I)*T+(G2(I)/T-G3(I))/T ELSE IF(I.EQ.4) THEN GAM=-95.23828+(62.656249-8.1454078*X)*X ELSE IF(I.EQ.5) THEN GAM=472.99219-74.144287*X-1869.6562/X2 ELSE IF(I.EQ.6) THEN GAM=825.17186-134.23096*X-2739.4375/X2 ELSE IF(I.EQ.7) THEN GAM=1181.3516-200.71191*X-2810.7812/X2 ELSE IF(I.EQ.8) THEN GAM=1440.1016-259.75781*X-1283.5625/X2 ELSE IF(I.EQ.9) THEN GAM=2492.1250-624.84375*X+30.101562*X2 ELSE IF(I.EQ.10) THEN GAM=4663.3129-1390.1250*X+97.671874*X2 ELSE GAM=I*I*I END IF COL(IT)=CT*EXP(-U0)*GAM ELSE IF(IC.GE.1) THEN GAM=I*I*I IF(I.LE.10) GAM=A(1,I)+A(2,I)*X+A(3,I)*X2+ * A(4,I)*X3+A(5,I)*X4+A(6,I)*X5 COL(IT)=CT*EXP(-U0)*GAM ELSE CALL CSPEC(II,NKI,IC,OSC0(IT),CPAR(IT),U0,T,COL(IT)) END IF C 100 I1=I+1 XI=I VI=XI*XI NHL=N1I-N0I+1 IF(N1Q.GT.0) NHL=N1Q IF(I1.GT.NHL) GO TO 200 C C ********** collisional excitation C C both explicit transitions as well as contributions to the C modified collisional ionization rate C DO 150 J=I1,NHL JJ=J+N0I-1 IC=0 IF(JJ.GT.N1I) GO TO 110 ICT=ITRA(II,JJ) IF(ICT.EQ.0) GO TO 150 IC=ICOL(ICT) 110 XJ=J VJ=XJ*XJ U0=ENION(N0I)*(1./VI-1./VJ)*TK IF(J.LE.20) C1=OSH(I,J) IF(J.GT.20) C1=OSH(I,20)*(20./XJ)**3 IF(IC.LT.0) GO TO 120 GAM=XI-(XI-1.)/(XJ-XI) IF(GAM.GT.XJ-XI) GAM=XJ-XI IF(I.GT.1) GAM=GAM*1.1 IF(U0.LE.UN) THEN EXPIU0=-LOG(U0)+EXPIA1+U0*(EXPIA2+U0*(EXPIA3+U0*(EXPIA4+ * U0*(EXPIA5+U0*EXPIA6)))) ELSE EXPIU0=EXP(-U0)*((EXPIB1+U0*(EXPIB2+U0*(EXPIB3+ * U0*(EXPIB4+U0))))/(EXPIC1+U0*(EXPIC2+ * U0*(EXPIC3+U0*(EXPIC4+U0)))))/U0 END IF CS=CT2/U0*C1*(0.693*EXP(-U0)+EXPIU0)*GAM GO TO 130 120 CALL CSPEC(II,JJ,IC,C1,CPAR(ICT),U0,T,COL(ICT)) GO TO 150 130 IF(JJ.GT.N1I) GO TO 140 COL(ICT)=CS GO TO 150 140 IF(IT.NE.0) COL(IT)=COL(IT)+CS 150 CONTINUE 200 CONTINUE RETURN END