OBNINSK_BCS.FOR


      PROGRAM BCSPHEN
      DIMENSION PAR(9),RES(7)
      ROBS(S,SIG,RO)=(2.*S+1.)/2./SIG*EXP(RO)
 1000 FORMAT(2X,2(F5.1,1X),8(F7.3,1X))
   10 READ 1000, Z,A,AAS,SHELL,D0,D1,OM2,BET,USTEP,UMAX
      PRINT 1000, Z,A,AAS,SHELL,D0,D1,OM2,BET,USTEP,UMAX
      PAR(1)=Z
      PAR(2)=A
      PAR(3)=AAS
C      PAR(3)=.073*A+.115*A**.66666
      PAR(4)=.40/A**.33333
      PAR(5)=SHELL
      PAR(6)=D0
      PAR(7)=D1
      PAR(8)=OM2
      PAR(9)=BET
C      PRINT 1000,(PAR(I),I=1,9)
      PRINT 1001
 1001 FORMAT(2X,' ENERGIES LOG.STATES & LEVELS    SIGMA2    TEMPER
     *QVIBR      QROT ')
 1010 FORMAT(2X,7(F8.3,2X))
 1020 FORMAT(2X,6(E10.4,1X))
      DU=USTEP
      NMAX=UMAX/USTEP
      DO 210 N=1,NMAX
      RES(1)=N*DU
      CALL BCS(PAR,RES)
      PRINT 1010,(RES(I),I=1,7)
  210 CONTINUE
      STOP
      END
      SUBROUTINE BCS(PAR,RES)
C*****PHENOMEN.DESCRIPTION BCS+COLL******
      DIMENSION PAR(9),RES(7)
      FU(X,DX,GX)=1.+(1.-EXP(-GX*X))*DX/X
      TM2=.24*PAR(2)**.66666
      Z=PAR(1)
      A=PAR(2)
      OM2=PAR(8)
      IF(OM2.LE.0.) OM2=30./A**.66666
      OM3=50./A**.66666
      BET=PAR(9)
      U=RES(1)
      IF(U-.01) 10,10,11
   10 DO 12 I=1,6
   12 RES(I+1)=0.
      RETURN
   11 AP=PAR(3)
      DL0=PAR(6)
      GAM=PAR(4)
      DW=PAR(5)
      NZ=Z
      NNZ=(NZ/2)*2
      IZ=0
      IF(NZ.GT.NNZ) IZ=1
      NN=A-Z
      NNN=(NN/2)*2
      IN=0
      IF(NN.GT.NNN) IN=1
      DEL=DL0*(IZ+IN)+PAR(7)
      U=U+DEL
      TKP=0.567*ABS(DL0)
      TKP2=TKP**2
      AX=AP*(1.+DW*GAM)
      DO 8 LX=1,5
      X=AX*TKP2
      AKP=AP*FU(X,DW,GAM)
      IF(ABS(AKP-AX)/AKP.LE.0.001) GO TO 9
      AX=AKP
    8 CONTINUE
    9 ECOH=0.152*AKP*DL0**2
      UKP=AKP*TKP2+ECOH
      SKP=2.*AKP*TKP
      DTKP=45.84*AKP**3*TKP**5
      FKP=0.608*AKP*TM2*(1.-0.6667*BET)
      FKP1=0.608*AKP*TM2*(1.+0.333*BET)
      IF(U.GE.UKP) GO TO 30
      IF(UKP.NE.0.) FI2=1-U/UKP
      FI=SQRT(FI2)
      T=2.*TKP*FI/ALOG((FI+1.)/(1.-FI))
      T1=T
      SK=SKP*TKP/T*(1-FI2)
      DET=DTKP*(1.-FI2)*(1.+FI2)**2
      FM=FKP*TKP/T*(1.-FI2)
      FM1=FKP1*(1.+2.*TKP/T*(1.-FI2))/3.
      GO TO 40
   30 UX=U-ECOH
      AS=AP*FU(UX,DW,GAM)
      IF(AS.LE.0.) RETURN
      T=SQRT(UX/AS)
      T1=T
      SK=2.*AS*T
      DET=45.84*AS**3*T1**5
      FM=0.608*AS*TM2*(1.-.6667*BET)
      FM1=0.608*AS*TM2*(1.+0.333*BET)
   40 SIG=FM*T1
      SIG1=FM1*T1
      IF(ABS(BET).GT.0.) SIG=FM1**.66666*FM**.33333*T1
      RES(2)=SK-ALOG(SQRT(DET))
      RES(3)=SK-ALOG(SQRT(6.283*SIG*DET))
      RES(4)=SIG
      RES(5)=T1
      CGA=.0075*A**.33333
      CALL QVIBR(T1,OM2,CGA,5,Q2)
      CALL QVIBR(T1,OM3,CGA,7,Q3)
      U=RES(1)
      CALL QROT(A,BET,SIG1,U,QR)
      RES(6)=Q2*Q3
      RES(7)=QR
      COL=ALOG(Q2*Q3*QR)
      RES(2)=RES(2)+COL
      RES(3)=RES(3)+COL
    2 RETURN
      END
C
      SUBROUTINE QVIBR(T,OM,CGA,LAM,Q)
C***** QVIBR INCLUDING DAMPING ***
      Q=1.
      IF(T.LT.0.05) GOTO 10
      GAM=CGA*(OM**2+(2.*3.141593*T)**2)
      FN=EXP(-GAM/OM/2)/(EXP(OM/T)-1.)
      IF(FN.LT.0.) GOTO 10
      U=LAM*OM*FN
      S=LAM*((1.+FN)*ALOG(1.+FN)-FN*ALOG(FN))
      Q=EXP(S-U/T)
   10 RETURN
      END
C
      SUBROUTINE QROT(A,BET,SIG,U,QR)
C***** QROT INCLUDING DAMPING ***
      FR(U)=1./(1.+EXP((U-UCR)/DCR))
      UCR=120.*BET*BET*A**.33333
      DCR=1400.*BET*BET/A**.66666
      IF(BET) 10,10,11
   10 QR=1
      GOTO 12
   11 QR=FR(U)*(SIG-1.)+1.
   12 RETURN
      END