SEMF code – Students assignment

program semf_assignment
implicit none

!—————————————————————————
! Declaring all the variables and constant involve in the SEMF equation
! Integer value input for A and Z since the number is always an integer.
!—————————————————————————
integer opt, A, Z
real a_vol, a_surf, a_coul, a_sym, a_pair, exp_pair
real p, q, r, s, t, BE
real BE_exp, BE_diff, Percentage

!—————————————————————————
! Providing interface to user for them to understand what is this program about
!—————————————————————————
write(*,*) “——————————————————————————————-”
write(*,*) “Calculation of Binding Energy Using Semi Empirical Mass Formula(SEMF) for EVEN-EVEN stable Nuclei”
write(*,*) “——————————————————————————————-”
print*,

!————————————————————————–
! Giving options to user so that they managed to choose the sets of
!coefficient they wanted to use in calculation
!————————————————————————–
write(*,*) “——————————————————————————————-”
print*, “Which sets of coefficient are you going to use in the calculation?”
print 101, “(1)”, “Least-Square fit 1”
print 101, “(2)”, “Least-Square fit 2”
print 101, “(3)”, “Rohlf”
print 101, “(4)”, “Wapstra”
print 101, “(5)”, “Own Coefficient”
101 format (a, x, a)
write(*,*) “——————————————————————————————-”

!—————————————————————————-
! Getting choices from user
!—————————————————————————-
read*, opt

!—————————————————————————-
! Asking user to input value of A and Z to perform SEMF calculation
!—————————————————————————-
print*, “What is the A value(mass number) of your nuclei?”
read*, A

print*, “What is the Z value(proton number) of your nuclei?”
read*, Z

!—————————————————————————-
! Sets of coefficient according to the selected option
!—————————————————————————-

if (opt == 1) then
a_vol = 15.8
a_surf = 18.3
a_coul = 0.714
a_sym = 23.2
a_pair = 12.0
exp_pair = 1.0/2.0

else if (opt == 2) then
a_vol = 15.76
a_surf = 17.81
a_coul = 0.711
a_sym = 23.702
a_pair = 34.0
exp_pair = 3.0/4.0

else if (opt == 3) then
a_vol = 15.75
a_surf = 17.8
a_coul = 0.711
a_sym = 23.7
a_pair = 11.18
exp_pair = 1.0/2.0

else if (opt == 4) then
a_vol = 14.1
a_surf = 13.0
a_coul = 0.595
a_sym = 19.0
a_pair = 33.5
exp_pair = 3.0/4.0

else if (opt == 5) then
print*, “What is the coefficient of the volume term?”
read*, a_vol

print*, “What is the coefficient of the surface term?”
read*, a_surf

print*, “What is the coefficient of the coulomb term?”
read*, a_coul

print*, “What is the coefficient of the asymmetry term?”
read*, a_sym

print*, “What is the coefficient of the pairing force?”
read*, a_pair

print*, “What is the power of A?”
read*, exp_pair

end if

!—————————————————————————-
! Equation for Semi Empirical Mass Formula
!—————————————————————————-

p = a_vol * A
q = a_surf * A **(2.0/3.0)
r = a_coul * Z * (Z-1.0)/ (A **(1.0/3.0))
s = a_sym * (A – (2.0* Z))**2.0 / A
t = a_pair / (A**(exp_pair))

BE = p – q – r – s + t

call Binding_Energy_Database(A, Z, BE_exp)

!—————————————————————————–
!The difference between the theoretical and experimental value is compute too
!—————————————————————————–

BE_diff = BE – (BE_exp * A)
Percentage = (BE_diff/ BE) * 100

!—————————————————————————–
! Print result to terminal and write to a file
!—————————————————————————–

print 201, “The volume term for this nuclei is “, p, “MeV”
print 201, “The surface term for this nuclei is “, q, “MeV”
print 201, “The coulomb term for this nuclei is “, r, “MeV”
print 201, “The asymmetry term for this nuclei is “, s, “MeV”
print 201, “The pairing force for this nuclei is “, t, “MeV”
print 201, “The binding energy for this nuclei is “, BE, “MeV”
print 201, “The difference between theoretical value and experimental data is “, BE_diff, “MeV”
print 201, “The percentage difference is “, Percentage, “%”
open (unit= 200, file= ‘result’, status= “replace”)
write (200,'(a, x, f9.4, x, a)’) “The binding energy for this nuclei is “, BE, “MeV”
write (200,'(a, x, f9.4, x, a)’) “The difference between theoretical value and experimental data is “, BE_diff, “MeV”
write (200,'(a, x, f9.4, x, a)’) “The percentage difference is “, Percentage, “%”

!—————————————————————————-
! Format of the Print Terminal
!—————————————————————————-
201 format(a, F9.4, a)

end program

!—————————————————————————-
! Subroutine for Database
!—————————————————————————-

subroutine Binding_Energy_Database(A, Z, BE_exp)
implicit none
integer A, Z
real, dimension(300,300):: BE_database
real BE_exp

BE_database(4,2) = 7.073915
BE_database(12,6) = 7.680144
BE_database(16,8) = 7.976206
BE_database(18,8) = 7.7670972
BE_database(20,8) = 7.568574
BE_database(20,10) = 8.03224
BE_database(22,10) = 8.080464810
BE_database(24,12) = 8.260709010
BE_database(26,12) = 8.333870110
BE_database(28,14) = 8.447744
BE_database(30,14) = 8.520654310
BE_database(32,16) = 8.493129
BE_database(34,16) = 8.583498010
BE_database(36,16) = 8.5753895
BE_database(36,18) = 8.519909210
BE_database(38,18) = 8.6142805
BE_database(40,18) = 8.595259
BE_database(42,20) = 8.6165634
BE_database(44,20) = 8.6581757
BE_database(46,20) = 8.668985
BE_database(46,22) = 8.6564514
BE_database(48,22) = 8.723005920
BE_database(50,22) = 8.755717820
BE_database(52,24) = 8.7759897
BE_database(54,24) = 8.7779557
BE_database(54,26) = 8.7363827
BE_database(56,26) = 8.7903545
BE_database(58,26) = 8.7922506
BE_database(58,28) = 8.7320596
BE_database(60,28) = 8.7807746
BE_database(62,28) = 8.7945537
BE_database(64,28) = 8.7774617
BE_database(66,30) = 8.75963211
BE_database(68,30) = 8.75568012
BE_database(70,32) = 8.72170012
BE_database(72,32) = 8.731745110
BE_database(74,32) = 8.7252
BE_database(76,32) = 8.705236
BE_database(74,34) = 8.687715
BE_database(76,34) = 8.711477
BE_database(78,34) = 8.717805720
BE_database(80,34) = 8.71081312
BE_database(82,34) = 8.6931966
BE_database(80,36) = 8.6929289
BE_database(82,36) = 8.710675
BE_database(84,36) = 8.717446
BE_database(84,38) = 8.67751215
BE_database(86,38) = 8.708456
BE_database(88,38) = 8.732595
BE_database(90,40) = 8.709968810
BE_database(92,40) = 8.692677710
BE_database(94,40) = 8.666800820
BE_database(92,42) = 8.657730520
BE_database(94,42) = 8.662333020
BE_database(96,42) = 8.653987310
BE_database(98,42) = 8.635168020
BE_database(100,44) = 8.6193593
BE_database(102,44) = 8.6074274
BE_database(104,44) = 8.58739924
BE_database(102,46) = 8.5802905
BE_database(104,46) = 8.58484813
BE_database(106,46) = 8.57999210
BE_database(108,46) = 8.56702310
BE_database(110,46) = 8.5471626
BE_database(110,48) = 8.5512753
BE_database(112,48) = 8.544730520
BE_database(114,50) = 8.522566
BE_database(116,50) = 8.523116210
BE_database(118,50) = 8.5165334
BE_database(120,50) = 8.5044927
BE_database(122,50) = 8.48790720
BE_database(120,52) = 8.477033
BE_database(122,52) = 8.47814012
BE_database(124,52) = 8.47327912
BE_database(126,52) = 8.46324812
BE_database(126,54) = 8.443543
BE_database(128,54) = 8.4432988
BE_database(130,54) = 8.437731
BE_database(132,54) = 8.427622
BE_database(130,56) = 8.40554920
BE_database(134,56) = 8.40817920
BE_database(136,56) = 8.402754920
BE_database(138,56) = 8.393419920
BE_database(140,58) = 8.37631711
BE_database(142,60) = 8.34603010
BE_database(146,60) = 8.3040929
BE_database(148,60) = 8.27717714
BE_database(144,62) = 8.30367911
BE_database(150,62) = 8.2616219
BE_database(152,62) = 8.2440618
BE_database(154,62) = 8.2268359
BE_database(154,64) = 8.2247968
BE_database(156,64) = 8.2153228
BE_database(158,64) = 8.2018198
BE_database(156,66) = 8.1924338
BE_database(158,66) = 8.19013015
BE_database(160,66) = 8.1840545
BE_database(162,66) = 8.1734575
BE_database(164,66) = 8.1587145
BE_database(162,68) = 8.1523975
BE_database(164,68) = 8.1490205
BE_database(166,68) = 8.1419597
BE_database(168,68) = 8.1296017
BE_database(170,68) = 8.1119599
BE_database(168,70) = 8.1118987
BE_database(170,70) = 8.106609
BE_database(172,70) = 8.097429
BE_database(174,70) = 8.083847
BE_database(176,70) = 8.064085
BE_database(176,72) = 8.0613598
BE_database(178,72) = 8.0494428
BE_database(180,72) = 8.0349308
BE_database(182,74) = 8.0183084
BE_database(184,74) = 8.0050774
BE_database(188,76) = 7.9738644
BE_database(190,76) = 7.9621043
BE_database(192,76) = 7.94852512
BE_database(192,78) = 7.94249113
BE_database(194,78) = 7.9359413
BE_database(196,78) = 7.9265293
BE_database(198,78) = 7.91415011
BE_database(196,80) = 7.91436915
BE_database(198,80) = 7.911551820
BE_database(200,80) = 7.9058953
BE_database(202,80) = 7.8968503
BE_database(204,80) = 7.885544920
BE_database(206,82) = 7.85753626
BE_database(208,82) = 7.8674536
BE_database(232,90) = 7.6150336
BE_database(238,92) = 7.5701256

BE_exp = BE_database(A,Z)

if (BE_exp == 0) then

print*, “Perhaps you may enter a nuclei which is not Even-Even or unstable nuclei.”
print*, “This application can only perform calculation for EVEN-EVEN stable nuclei.”
print*, “There is no experimental value for this nuclei.”
print*, “Try it again with a more suitable nuclei. GoodLuck!!”

stop

end if

end subroutine