MODULE functs USE gem_kind USE trig IMPLICIT NONE PRIVATE PUBLIC :: gauss, cl_go, factorial, exp_decay, transform, ran_j PUBLIC :: get_little_a, OPERATOR(+), zbeta PRIVATE :: addvectors INTERFACE OPERATOR (+) MODULE PROCEDURE addvectors END INTERFACE REAL (KIND=r8), DIMENSION(0:300), PUBLIC :: FCT REAL (KIND=r8), PARAMETER, PRIVATE :: fact=50.0 CONTAINS FUNCTION gauss () RESULT (FUNCT_out) !+ ! ! FUNCTIONAL DESCRIPTION: ! ! gaussian random number generator centered ! on zeroed with standard deviation ! of unity ! !- REAL (kind=r4) :: funct_out REAL (kind=r4), SAVE :: x=0.0,theta=0.0 REAL (KIND=r4) :: ran LOGICAL, SAVE :: first = .true. IF (first) THEN CALL RANDOM_NUMBER(ran) x = (-2.0*LOG(RAN+1.0e-37))**0.5 CALL RANDOM_NUMBER(ran) theta = 360.0*RAN funct_out = x*COSD(theta) first = .FALSE. ELSE funct_out = x*SIND(theta) first = .TRUE. END IF END FUNCTION gauss PURE FUNCTION cl_go(J1,J2,J,M1,M2,M) RESULT(funct_out) !calculates clebsch Gordan coeeficients ! J's and M's must be twice the real values INTEGER, INTENT(IN) :: J1,J2,J,M1,M2,M REAL (kind=r4) :: funct_out INTEGER :: Ja,Jb,Jc,Jd,Je,Jf,Jg,Jh,Ji,Jj INTEGER :: fase,zmin,zmax,z REAL (KIND=r8) :: cg,f1,f2 funct_out = 0.0 IF (M1+M2 /= M) RETURN IF (ABS(M1) > J1) RETURN IF (ABS(M2) > J2) RETURN IF (ABS(M) > J) RETURN IF (J > J1+J2) RETURN IF (J < ABS(J1-J2)) RETURN IF (j == 0 .AND. J1 == 0 .AND. j2 == 0)THEN funct_out = 1.0 RETURN END IF Zmin = MAX(0, -J+J2-M1, -J+J1+M2) Zmax = MIN(J1+J2-J, J2+M2, J1-M1) cg = 0.0 DO z=zmin,zmax,2 Ja = Z/2 Jb = (J1 + J2 - J - Z)/2 Jc = (J1 - M1 - Z)/2 Jd = (J2 + M2 - Z)/2 Je = (J - J2 + M1 + Z)/2 Jf = (J - J1 - M2 + Z)/2 fase = ((-1)**(z/2)) f2 = REAL(fase, KIND=r8) cg = cg + f2/(FCT(Ja)*FCT(Jb)*FCT(Jc)*FCT(Jd)*FCT(Je)*FCT(Jf)) END DO Ja = (J1 + J2 - J)/2 Jb = (J1 - J2 + J)/2 Jc = (-J1 + J2 + J)/2 Jd = (J1 + M1)/2 Je = (J1 - M1)/2 Jf = (J2 + M2)/2 Jg = (J2 - M2)/2 Jh = (J + M)/2 Ji = (J - M)/2 Jj = (J1 + J2 + J + 2)/2 F1 = REAL(J + 1, KIND=r8) cg = SQRT(F1*FCT(Ja)*FCT(Jb)*FCT(Jc)*FCT(Jd)*FCT(Je)*FCT(Jf)* & FCT(Jg)*FCT(Jh)*FCT(Ji)/FCT(Jj))*cg funct_out = REAL(cg/fact**0.5) RETURN END FUNCTION cl_go SUBROUTINE factorial () INTEGER :: i FCT(0) = 1.0_r8 FCT(1) = 1.0_r8/fact DO i=2,300 FCT(i) = FCT(i-1)*REAL(i, KIND=r8)/fact END DO RETURN END SUBROUTINE factorial FUNCTION exp_decay(decay_width) RESULT(funct_out) ! for a decay width in MEV, returns in zs, ! a random time appropriate for exponential decay ! law REAL (kind=r4), INTENT(IN) :: decay_width REAL (kind=r4) :: funct_out REAL (kind=r4) :: ran CALL RANDOM_NUMBER(ran) funct_out = - 0.65824*LOG(RAN+1.0e-37)/decay_width RETURN END FUNCTION exp_decay PURE SUBROUTINE transform(theta1,phi1,theta2,phi2,theta3,phi3) ! ! FUNCTIONAL DESCRIPTION: ! ! this subroutine performs a rotational transformation. ! theta1,phi1 are the coordinates (spherical angles in ! degrees)of a unit vector in the original coordinate systems. ! the z axis is made to rotate in the phi=phi2 plane by an ! angle theta2. The coordinates of the vector in the new ! reference frame are theta3,phi3 ! ! ! To calculate the angle between two vectors ! theta1,phi1 and theta2,phi2 ! use CALL transform(theta1,phi1,-theta2,phi2,theta3,phi3) ! i.e. theta3 in the angle between them, note the negative sign ! on theta2 ! ! ! ! DUMMY ARGUMENTS: ! ! input: theta1,phi1 (REAL (KIND=r4) ::) theta and phi angles (degrees) of initial unit ! vector ! ! theta2,phi2 (REAL (KIND=r4) ::) rotational angles ! ! OUTPUT: theta3,phi3 (REAL (KIND=r4) ::) theta and phi angles(degrees) of unit vector ! after rotational transform ! REAL (KIND=r4), INTENT(IN) :: theta2,phi2 REAL (KIND=r4), INTENT(IN) :: theta1,phi1 REAL (KIND=r4), INTENT(OUT) :: theta3,phi3 REAL (KIND=r4) :: xp,yp,zp REAL (KIND=r4) :: xt,yt,zt ! rotate vector in x-y plane by -phi2 ! and find cartesian coordinates xp = SIND(theta1)*COSD(phi1-phi2) yp = SIND(theta1)*SIND(phi1-phi2) zp = COSD(theta1) ! rotate vector in z-x plane by theta2 zt = zp*COSD(theta2) - xp*SIND(theta2) xt = xp*COSD(theta2) + zp*SIND(theta2) yt = yp ! rotate vector in x-y plane back by +phi2 ! and find spherical coordinates IF (xt == 0.0 .AND. yt == 0.0) THEN phi3 = 0 ELSE phi3 = ATAN2d(yt,xt) + phi2 if (phi3 > 360.0) phi3 = phi3 - 360.0 if (phi3 < 0.0) phi3 = phi3 + 360.0 END IF if (zt > 1.0) THEN theta3 = 0 ELSE IF(zt < -1.0) THEN theta3 = 180 ELSE theta3 = ACOSD(zt) END IF RETURN END SUBROUTINE transform PURE FUNCTION addvectors (v1,v2) RESULT(v3) ! ! FUNCTIONAL DESCRIPTION: ! ! The subroutine adds together two vectors "1"+"2" specified in spherical ! polar coordinates also giving the answer "3" in polar coordinates ! ! DUMMY ARGUMENTS: ! ! input: r1,theta1,phi1 (type (vell)) magnitude, theta and phi angle (degrees) of ! vector "1" ! r2,theta2,phi2 (type (vell)) - vector "2" ! ! output: r3,theta3,phi3 (type (vell) vector "3" = "1" + "2" TYPE (vell), INTENT(IN) :: v1,v2 TYPE (vell) :: v3 REAL (KIND=r4) :: x,y,z z = v1%v*COSD(v1%theta) + v2%v*COSD(v2%theta) x = v1%v*SIND(v1%theta)*COSD(v1%phi) + v2%v*SIND(v2%theta)*COSD(v2%phi) y = v1%v*SIND(v1%theta)*SIND(v1%phi) + v2%v*SIND(v2%theta)*SIND(v2%phi) v3%v = (x**2.0+y**2.0+z**2.0)**0.5 IF (v3%v == 0.0) THEN v3%theta = 0.0 v3%phi = 0.0 ELSE v3%theta = ACOSD(z/v3%v) IF (y == 0.0 .AND. x == 0.0) THEN v3%phi = 0.0 ELSE v3%phi = ATAN2D(y,x) IF (v3%phi < 0.0) v3%phi = v3%phi + 360.0 END IF END IF RETURN END FUNCTION addvectors SUBROUTINE ran_j(x,ac,U,aden,rua,xmax) ! ! FUNCTIONAL DESCRIPTION: ! ! this subroutine chooses in the Monte Carlo fashion the angular momentum ! associated with one or other of the angular momentum bearing modes; i.e. ! twisting, tilting, bending etc. ! ! DUMMY ARGUMENTS: ! ! input: ac (REAL (KIND=r4) ::) energy of modes = ac*x**2./2. MeV ! U (REAL (KIND=r4) ::) - thermal excitation energy at saddle-point ! aden (REAL (KIND=r4) ::) - little "a" level density parameter MeV-1 ! rua (REAL (KIND=r4) ::) - 2.*(aden*U)**.5 ! xmax (REAL (KIND=r4) ::) - maximum value of x possible, i.e. for the tilting ! mode the maximum projection is the total angular momentum ! ! output: x angular momentum (hbar) ! ! IMPLICIT INPUTS: ! ! none ! ! IMPLICIT OUTPUTS: ! ! none ! ! ! SIDE EFFECTS: ! ! none ! ! REAL (KIND=r4) :: y1,y2,xmax1,xmax2,ran REAL (KIND=r4), INTENT(OUT) :: x REAL (KIND=r4), INTENT(IN) :: aden,rua,xmax,U,ac xmax1 = (2.0*U/ac)**0.5 xmax2 = MIN(xmax,xmax1) DO CALL RANDOM_NUMBER(ran) x = xmax2*(1.0-2.0*RAN) y1 = exp(2.0*(aden*(U-0.5*ac*x**2.0))**0.5-rua) CALL RANDOM_NUMBER(y2) IF (y2 < y1) EXIT END DO RETURN END subroutine ran_j SUBROUTINE get_little_a(A,Z,J,U,surf,shell_corr,E_rot,aden,entropy,& pair_energy,aa_scale) !+ ! ! calculates the level density parameter "litte a" aden ! aden = little a ! aden_1 = little at zero temperature, if little a is ! temperature dependent ! !- REAL (KIND=r4), INTENT(IN) :: A,Z,J,surf,aa_scale,pair_energy,shell_corr,E_rot REAL (KIND=r4), INTENT(INOUT) :: U REAL (KIND=r4), INTENT(OUT) :: aden,entropy REAL (KIND=r4) :: aden_1,aden_2,us,facte,E_sym,kappa REAL (KIND=r4), PARAMETER :: E_fade = 18.52 INTEGER :: i !backshift excitation energy U = U + pair_energy IF (para%mass_option == 0) THEN U = U + shell_corr !liquid drop mass is used ELSE IF (para%mass_option == 2) THEN u = u + shell_corr*(1.0-EXP(-(U+E_rot)/e_fade)) !shell corrections !fade out ! u = u + shell_corr/(1.0+EXP(-(U+E_rot-55.)/5.)) !shell corrections !fade out ! u = u + shell_corr*(1.0-EXP(-U/e_fade)) !shell corrections !fade out END IF E_sym = 0.0 SELECT CASE (para%aden_type) CASE (0) aden = A/para%Aden_0 CASE (-1) !Toke + Swiatecki - one sphere aden = A/14.61*(1+4.0*surf/A**0.3333333) CASE (-2) !Ignatyuk et al aden = A*0.073 + 0.095*surf*A**0.6666666 ! aden = A*0.073 + 0.095*A**0.6666666 CASE (-3) !Gottschalk + Ledergerber aden = A*0.0999373 - 0.045934*surf*A**0.6666666 CASE (-4) ! Ormand et al temp dependent aden_1 = A/8.0 aden = aden_1/1.4*(1.0+0.4*EXP(0.1-1.213*U/A)) case (-5) aden_1 = A/8.0& *((Z/A)**0.33333333+((A-z)/A)**0.333333333)/2.0**0.3333333333 aden = aden_1/1.4*(1.0+0.4*EXP(0.1-1.213*U/A)) CASE (-6) ! lestone ! aden_1 = (1.6*A+1.8*surf*A**0.666666666)/15.5 aden_1 = (1.6*A+1.8*A**0.666666666)/15.5 aden_2 = 0.5*A/15.5 ! no temperature dependent symmerty energy ! us = (U-E_sym)/441.0/A**0.33333333 ! facte = 1.0-EXP(-Us*8.0) ! us = Us*A/(aden_1-aden_2*facte) ! aden = aden_1 - aden_2*(1.0-EXP(-us)) ! with temperature dependent symmetry energy DO i=1,2 us = (U-E_sym)/441/A**0.33333333 facte = 1.0-EXP(-Us*8.0) us = Us*A/(aden_1-aden_2*facte) aden = aden_1 - aden_2*(1.0-EXP(-us)) E_sym = 0.82247*(1./aden-1.0/aden_1)*(A-2*Z)**2.0 END DO U = U - E_sym CASE (-7) aden = A/11.5 + 0.117*surf*A**0.66666*EXP(-U/A*5/3.0) CASE (-8) aden = A/15 + 0.34*surf*A**0.66666*EXP(-U/A*5/3.0) CASE (-9) !Fineman aden = A/(8.0+para%aden_0*U/A) CASE (-10) aden = A/(7.0 + para%aden_0*U/A) CASE (-11) aden = A/(6.0 + para%aden_0*U/A) CASE (-12) aden = A/(8.0+para%aden_0*(U/A)**2.) CASE (-13) aden = (a-z)/para%aden_0 CASE (-14) aden = z/para%aden_0 CASE (-15) !grimes CASE C aden = A/8.7873*EXP(-0.0493*(z-zbeta(A))**2.) CASE (-16) !grimes case B aden = A/9.009*EXP(-0.000641*(A-2*z)**2.) CASE (-17) !grimes case B aden = A/12.024*EXP(0.000461*(A-2*z)**2.) CASE (-18) !grimes case B ! aden = A/8*(1.-0.1*U/A+0.01*(A-2*z)*U/A) ! aden = A/8*(1.-0.2*U/A+0.01*(A-2*z)*U/A) ! aden = A/8*(1.-0.3*U/A+0.02*(A-2*z)*U/A) ! aden = A/8*(1.-0.4*U/A+0.02*(A-2*z)*U/A) ! aden = A/8.*(1-0.0125*(A-2*z)*exp(-u/A/0.31)) ! aden = A/10.*(1-0.02*(A-2*z)*exp(-u/A/0.31)) ! aden = A/12.*(1-0.02*(A-2*z)*exp(-u/A/0.31)) ! aden = A/12.-A/12.*0.03*(A-2*z)*exp(-u/A/0.31)) ! aden = A/12.-A/12.*0.04*(A-2*z)*exp(-u/A/0.31) ! aden = A/12.-A/12.*0.04*(A-2*z)*exp(-u/A/0.51) ! aden = A/12.-A/12.*0.04*(A-2*z)*exp(-u/A/0.21) ! aden = A/12. -0.5*(A-2.*Z)*exp(-u/A/0.31) + 0.2*(A-2.*Z) ! aden = A/14. -0.5*(A-2.*Z)*exp(-u/A/0.31) + 0.2*(A-2.*Z) ! aden = A/14. -0.5*(A-2.*Z)*exp(-u/A/0.51) + 0.2*(A-2.*Z) ! aden = A/16. -0.5*(A-2.*Z)*exp(-u/A/0.51) + 0.2*(A-2.*Z) ! aden = A/16. -0.4*(A-2.*Z)*exp(-u/A/0.81) + 0.2*(A-2.*Z) ! good ! aden = 0.075*A + 0.0625*A*EXP(-U/A/0.51) & ! -0.5*(A-2.*Z)*exp(-u/A/0.51) + 0.2*(A-2.*Z) ! aden = 0.075*A + 0.075*A*EXP(-U/A/0.51) & ! -0.6*(A-2.*Z)*exp(-u/A/0.51) + 0.2*(A-2.*Z) !fig25 ! aden = 0.069*A + 0.081*A*EXP(-U/A/0.6) & ! -0.65*(A-2.*Z)*exp(-u/A/0.6) + 0.25*(A-2.*Z) ! aden = 0.1005*A + 0.081*A*EXP(-U/A/0.4) & ! -0.65*(A-2.*Z)*exp(-u/A/0.4) !fig27 ! kappa = 0.2*(Z-zbeta(A)) ! aden = A/(8.+ kappa*U/A) !fig28 kappa = 0.1*(z-zbeta(A)) ! kappa = 0.1*(z-zbeta(A))-.4 !fig29 ! kappa= 0.4 ! kappa = 8.2e-5*Z**2 ! kappa = 1.6e-4*z**2 ! kappa = 0.6 + (z-70.)*0.1 aden = A/(8.+ kappa*(U/A)**2.) CASE (-19) !grimes case B aden = A**(2./3.)*(Z**(1./3)+(A-Z)**(1./3.))/2.**(2./3.)& /para%aden_0 CASE (-20) !my calculations aden = 0.861 + 0.08361*A + (A-2.*Z)**2*(6.229e-2/A-8.23e-5) CASE (-21) !my calculations aden = A/(8.+para%aden_0*U/A*J/50.) END SELECT aden = aden*aa_scale IF (U > 0.0) THEN entropy = 2.0*(Aden*U)**0.5 ELSE entropy = 0.0 END IF RETURN END SUBROUTINE get_little_a ELEMENTAL FUNCTION zbeta(A) RESULT(funct_out) REAL (KIND=r4), INTENT(IN) :: A REAL (KIND=r4) :: funct_out funct_out = A/(1.98+0.0155*A**(2./3.)) RETURN END FUNCTION zbeta END MODULE functs