IMPLICIT REAL*8(a-h,o-z) CHARACTER*29 FRASA CHARACTER*23 FRASE PARAMETER(nbox=20, zero=0.0D0) DIMENSION iopab(nbox),iopbc(nbox),iopac(nbox), iopsum(nbox), . iopsumtot(nbox) DIMENSION opab(nbox),opbc(nbox),opac(nbox), opsum(nbox), . opsumtot(nbox) DIMENSION idangab(nbox),idangbc(nbox),idangac(nbox), . idangsum(nbox) DIMENSION dangab(nbox),dangbc(nbox),dangac(nbox),dangsum(nbox) DIMENSION ietrab(nbox),ietrbc(nbox),ietrac(nbox), ietrsum(nbox) DIMENSION etrab(nbox), etrbc(nbox), etrac(nbox), etrsum(nbox) DIMENSION ievibc(nbox), ieviab(nbox), ieviac(nbox), ievisum(nbox) DIMENSION evibc(nbox), eviab(nbox), eviac(nbox), evisum(nbox) DIMENSION ierobc(nbox), ieroab(nbox), ieroac(nbox), ierosum(nbox) DIMENSION erobc(nbox), eroab(nbox), eroac(nbox), erosum(nbox) DIMENSION angle(nbox), anglebin(nbox+1), pimp (nbox+1), . etab(nbox+1), etbc(nbox+1), etac(nbox+1), . evab(nbox), evbc(nbox), evac(nbox), . erab(nbox), erbc(nbox), erac(nbox) PARAMETER (pi=3.141592654D0) OPEN(UNIT=10, FILE="input-data", STATUS="UNKNOWN") OPEN(UNIT=22, FILE="OPA-TOT", STATUS="UNKNOWN", ACCESS="APPEND") OPEN(UNIT=23, FILE="DC-TOT", STATUS="UNKNOWN", ACCESS="APPEND") OPEN(UNIT=24, FILE="ENER-TOT", STATUS="UNKNOWN", ACCESS="APPEND") WRITE(*,*)"Enter the total number of Tj" READ(*,*)ntj READ(10,*) READ(10,10)FRASB, bmax 10 FORMAT (A29,F10.5) READ(10,11)FRASE, etr 11 FORMAT (A23,F12.5) jump=6 DO i=1,jump READ(10,*) END DO C:: Read all the tjs. to calculate average energy distribution evib_sum= 0.0D0 erot_sum = 0.0D0 etrans_sum = 0.0D0 icount=0 DO i=1,ntj !ntj READ(10,101) k, ic, theta, evib, erot, eint, etrans, . epot, sang, ener_diff IF(ic.EQ.1.OR.ic.EQ.2) THEN evib_sum = evib_sum + evib erot_sum = erot_sum + erot etrans_sum = etrans_sum + etrans icount = icount +1 END IF END DO evib_aver = evib_sum/real(icount) erot_aver = erot_sum/real(icount) etrans_aver = etrans_sum/real(icount) write(98,*) evib_aver, erot_aver, etrans_aver 101 FORMAT(I5,9X,I1,6X,6(F9.4,1X),F9.4,4X,F9.4) C---------- C STOP C---------- C:: inizialization of arrays DO i=1,20 iopab(i)=0 iopbc(i)=0 iopac(i)=0 iopsum(i)=0 END DO C:: C:: reading of channel A+BC C:: DO i=1,6 READ(10,*) !jump END DO READ(10,180)emaxbc DO i=1,2 READ(10,*) !jump END DO DO l=1,20 read(10,190) iopbc(l), a1, ietrbc(l),ievibc(l),ierobc(l), . idangbc(l) write(99,*) l, ierobc(l) END DO DO i=1,4 READ(10,*) !jump END DO C---------- C STOP C---------- C:: C:: reading of channel C+AB C:: READ(10,180)emaxab DO i=1,2 READ(10,*) !jump END DO DO l=1,20 read(10,190) iopab(l), a1, ietrab(l),ieviab(l),ieroab(l), . idangab(l) END DO DO i=1,4 READ(10,*) END DO C---------- C STOP C---------- C:: C:: reading of channel B+AC C:: READ(10,180)emaxac DO i=1,2 READ(10,*) !jump END DO DO l=1,20 read(10,190) iopac(l), a1, ietrac(l), ieviac(l), ieroac(l), . idangac(l) END DO C---------- C STOP C---------- C:: C:: ranges for the distributions C:: C:: angle distribution DO k=1,20 angle(k)=real(k-1)*9.0+9.0/2.0 END DO C:: impact parameter distribution pimp(1) = 0.0D0 DO k=2,21 pimp(k)= real(k-2)*(bmax/20.0D0) + bmax/20.0D0 END DO C:: scattering angle distribution anglebin(1) = 0.0D0 DO k=2,21 anglebin(k)= real(k-2)*(180.0D0/20.0D0) + 180.0d0/20.0d0 END DO C:: energy distribution etab(1) = 0.0D0 etbc(1) = 0.0D0 etac(1) = 0.0D0 DO k=2,21 etab(k)=real(k-2)*emaxab/20.0+emaxab/20.0 etbc(k)=real(k-2)*emaxbc/20.0+emaxbc/20.0 etac(k)=real(k-2)*emaxac/20.0+emaxac/20.0 END DO C:: C:: Normalization of translational distribution C:: DO k=1,20 etrab(k)=real(ietrab(k)) etrbc(k)=real(ietrbc(k)) etrac(k)=real(ietrac(k)) etrsum(k)=etrab(k)/ntj + etrac(k)/ntj END DO etrabm=etrab(1) etrbcm=etrbc(1) etracm=etrac(1) etrsumm=etrsum(1) DO k=1,20 IF(etrab(k).GT.etrabm) etrabm=etrab(k) IF(etrbc(k).GT.etrbcm) etrbcm=etrbc(k) IF(etrac(k).GT.etracm) etracm=etrac(k) IF(etrsum(k).GT.etrsumm) etrsumm=etrsum(k) END DO DO k=1,20 IF(etrabm.NE.0.0) then etrab(k)=etrab(k)/etrabm etrbc(k)=etrbc(k)/etrbcm etrac(k)=etrac(k)/etracm etrsum(k)=etrsum(k)/etrsumm ENDIF END DO C:: Normalization of vibrational distribution DO k=1,20 eviab(k)=real(ieviab(k)) evibc(k)=real(ievibc(k)) eviac(k)=real(ieviac(k)) evisum(k)=eviab(k)/ntj + eviac(k)/ntj END DO eviabm=eviab(1) evibcm=evibc(1) eviacm=eviac(1) evisumm=evisum(1) DO k=1,20 IF(eviab(k).GT.eviabm) eviabm=eviab(k) IF(evibc(k).GT.evibcm) evibcm=evibc(k) IF(eviac(k).GT.eviacm) eviacm=eviac(k) IF(evisum(k).GT.evisumm) evisumm=evisum(k) END DO DO k=1,20 IF(etrabm.NE.0.0) then eviab(k)=eviab(k)/eviabm evibc(k)=evibc(k)/evibcm eviac(k)=eviac(k)/eviacm evisum(k)=evisum(k)/evisumm ENDIF END DO C:: Normalization of rotational distribution DO k=1,20 eroab(k)=real(ieroab(k)) erobc(k)=real(ierobc(k)) eroac(k)=real(ieroac(k)) erosum(k)=eroab(k)/ntj + eroac(k)/ntj END DO eroabm=eroab(1) erobcm=erobc(1) eroacm=eroac(1) erosumm=erosum(1) DO k=1,20 IF(eroab(k).GT.eroabm) eroabm=eroab(k) IF(erobc(k).GT.erobcm) erobcm=erobc(k) IF(eroac(k).GT.eroacm) eroacm=eroac(k) IF(erosum(k).GT.erosumm) erosumm=erosum(k) END DO DO k=1,20 IF(etrabm.NE.0.0) then eroab(k)=eroab(k)/eroabm erobc(k)=erobc(k)/erobcm eroac(k)=eroac(k)/eroacm erosum(k)=erosum(k)/erosumm ENDIF END DO C:: C:: Calculation of differential cross section C:: DO k=1,20 idangsum(k)=idangab(k) + idangac(k) END DO DO k=1,20 dangab(k)=real(idangab(k)) dangbc(k)=real(idangbc(k)) dangac(k)=real(idangac(k)) dangsum(k)=real(idangsum(k)) END DO DO k=1,20 dangab(k)=dangab(k)/sin(angle(k)*pi/180.0) dangbc(k)=dangbc(k)/sin(angle(k)*pi/180.0) dangac(k)=dangac(k)/sin(angle(k)*pi/180.0) dangsum(k)=pi*bmax*bmax*dangsum(k) / . (ntj*2.0D0*pi*(sin(angle(k)*pi/180.0))) END DO dangabm=dangab(1) dangbcm=dangbc(1) dangacm=dangac(1) DO k=1,20 IF(dangab(k).GT.dangabm) dangabm=dangab(k) IF(dangbc(k).GT.dangbcm) dangbcm=dangbc(k) IF(dangac(k).GT.dangacm) dangacm=dangac(k) END DO DO k=1,20 dangab(k)=dangab(k)/dangabm dangbc(k)=dangbc(k)/dangbcm dangac(k)=dangac(k)/dangacm END DO C:: sum the opacity function for the two equivalent channels DO l=1,20 iopsum(l) = iopab(l) + iopac(l) END DO C:: C:: Normalization of opacity function C:: DO l=1,20 iopsum(l) = iopab(l) + iopac(l) iopsumtot(l) = iopab(l) + iopac(l) + iopbc(l) opab(l)=real(iopab(l)) opbc(l)=real(iopbc(l)) opac(l)=real(iopac(l)) opsum(l)=real(iopsum(l)) opsumtot(l)=real(iopsumtot(l)) END DO opabm=opab(1) opbcm=opbc(1) opacm=opac(1) opsumm=opsum(1) DO k=1,20 IF(opab(k).gt.opabm) opabm=opab(k) IF(opbc(k).gt.opbcm) opbcm=opbc(k) IF(opac(k).gt.opacm) opacm=opac(k) IF(opsum(k).gt.opsumm) opsumm=opsum(k) END DO DO k=1,20 opab(k)=opab(k)/opabm opbc(k)=opbc(k)/opbcm opac(k)=opab(k)/opacm opsum(k)=opsum(k)/opsumtot(k) END DO C:: C:: Writing the histograms for channel SUM (AB and AC channels) C:: WRITE(22, 208) etr WRITE(22, 209) WRITE(23, 208) etr WRITE(23, 209) WRITE(24, 208) etr WRITE(24, 209) C:: Write Opacity functions WRITE(22, 204) DO k=1,20 write(22,201)pimp(k),opsum(k) END DO write(22,201)pimp(21),zero C:: Write differential cross section WRITE(23,210) DO k=1,20 write(23,201) anglebin(k),dangsum(k) END DO WRITE(23,201) anglebin(21),zero C:: Write energy distribution and averages C:: Write Etr WRITE(24,211) DO k=1,20 write(24,201)etac(k),etrsum(k) END DO WRITE(24,201) etac(21),zero WRITE(24,*)'---- TRANS. ENERGY AVERAGE= ', etrans_aver, ' ----' C:: Write Evib WRITE(24,212) DO k=1,20 write(24,201)etac(k),evisum(k) END DO WRITE(24,201) etac(21),zero WRITE(24,*)'---- EVIB. ENERGY AVERAGE= ', evib_aver, ' ----' C:: Write Erot WRITE(24,213) DO k=1,20 write(24,201)etac(k),erosum(k) END DO WRITE(24,201) etac(21),zero WRITE(24,*)'---- EROT. ENERGY AVERAGE= ', erot_aver, ' ----' 180 FORMAT(50x,F9.4) C 190 FORMAT(3x,4x,I5,8x,I5,10x,I5) 190 FORMAT(3x,4x,I5,4X,5(4X,I5,6X)) 199 FORMAT(/' ************** CHANNEL A+BC ************* ',/) 200 FORMAT(1x,'#TRANS. ENERGY DISTRIBUTION (normalizad)',/) 300 FORMAT(/,1x,'#VIB. ENERGY DISTRIBUTION (normalizad)',/) 400 FORMAT(/,1x,'#ROT. ENERGY DISTRIBUTION (normalizad)',/) 201 FORMAT(1x,f7.2,3x,f10.4) 202 FORMAT(/,1x,'#ANGULAR DISTRIBUTION (normalizad)',/) 203 FORMAT(1x,f7.2,3x,f5.3,/) 204 FORMAT(/,1x,'#OPACITY FUNCTION (normalized)',/) 205 FORMAT(1x,f7.2,3x,f5.3) 206 FORMAT(/,' ************** CHANNEL C+AB ************* ',/) 207 FORMAT(/,' ************** CHANNEL B+AC ************* ',/) 208 FORMAT(/,' Results for E_tr=', F12.5) 209 FORMAT(/,' ************** CHANNEL SUM (B+AC plus C+AB) . ************* ') 210 FORMAT(/,' DIFF. CROSS SECTION' ,/) 211 FORMAT(/,' TRANS. ENERGY DISTRIBUTION' ,/) 212 FORMAT(/,' VIB. ENERGY DISTRIBUTION' ,/) 213 FORMAT(/,' ROT. ENERGY DISTRIBUTION' ,/) END