!This program generates NE (set below) EoS candidates. The number of !points (the grid) in the energy-density/pressure plane can be set !below (Npr, Npe). The program outputs three files. The first one ! is called boundary_pts.dat. This file contains the coordinates of !the boundaries dictated by causality and monotony for a given choice !of boundary conditions. The number of points used for the boundary !is set by the parameter Nbd. The second output file is called !EoS_grid.dat. It contains the e-density and pressure coordinates !of the grid that is used to generate the EoS candidates. Finally, !a file EoS_points.dat is created. It has NE+1 lines. The first line !contains the e-ensity coordinates. The second line corresponds to !the coordinates of the first EoS candidates, the third line !to the second EoS candidate, and so on. Each line has !Ned+2 entries. The extra '2' comes from the start and end point !of a candidate, that is, the boundary value, which is chosen at !random (for now, only the pressure value is picked at random). PROGRAM main USE aux IMPLICIT NONE integer, parameter :: Nd=20, Np=20 !in low band, number of pts en densty (Nd),number of pts pressure (Np) integer, parameter :: Nep=200 !number of EoS integer, parameter :: Neh=12!number of points of density and pressure in high (pQCD) band double precision :: h1,E1,P0,Emax,Pmin,Pmax!sp1 double precision :: h,E(Nd),P(Nd), edensity (Nd)!Bandmax(Nd,2),Bandmin(Nd,2) double precision :: Pin(Nep)!,sp2(Nep) double precision :: Bandground (Nd),Bandroof (Nd),Bandfloor (Nd) double precision :: table(Neh,3), Bandhighground (Neh),Bandhighroof (Neh),Bandhighfloor(Neh), Ehigh(Neh),Phigh(Neh) double precision :: factor INTEGER(KIND=4), PARAMETER :: NE = 1 !number of EoS candidates to be generated INTEGER(KIND=4), PARAMETER :: Npr = 20 !number of pts pressure INTEGER(KIND=4), PARAMETER :: Ned = 20 !number of pts en densty INTEGER(KIND=4), PARAMETER :: Nbd = 1000 !number of pts for boundary plot REAL(KIND=8) :: P1x,P2x,P1y,P2y,P1intx, P1inty,P2intx,P2inty REAL(KIND=8), PARAMETER :: Q1x=15231.618285 !bndry pt (energy density at baryon chemical potential 2.6 GeV) REAL(KIND=8), PARAMETER :: Q1y=3465.932101 !bndry pt (lower pressure at baryon chemical ppotential 2.6 GeV) REAL(KIND=8), PARAMETER :: Q2x=15231.618285 !bndry pt (energy density at baryon chemical potential 2.6 GeV) REAL(KIND=8), PARAMETER :: Q2y=4642.5168 !bndry pt (higher pressure at baryon chemical potential 2.6 GeV) ! REAL(KIND=8), PARAMETER :: Q1x=20305.13 !bndry pt (energy density at baryon chemical potential 2.8 GeV) ! REAL(KIND=8), PARAMETER :: Q1y= 5063.841 !bndry pt (lower pressure at baryon chemical ppotential 2.8 GeV) ! REAL(KIND=8), PARAMETER :: Q2x=20305.13 !bndry pt (energy density at baryon chemical potential 2.8 GeV) ! REAL(KIND=8), PARAMETER :: Q2y=6244.83 !bndry pt (higher pressure at baryon chemical potential 2.8 GeV) ! REAL(KIND=8), PARAMETER :: P1intx=P1y-Q1y+Q1x !bndry pt ! REAL(KIND=8), PARAMETER :: P1inty=P1y !bndry pt ! REAL(KIND=8), PARAMETER :: P2intx=Q2y-P2y+P2x !bndry pt ! REAL(KIND=8), PARAMETER :: P2inty=Q2y !bndry pt ! REAL(KIND=8) :: Q1x,Q1y,Q2x,Q2y REAL(KIND=8) :: P1(2),P2(2),Q1(2),Q2(2),P1int(2),P2int(2) !bndry pts REAL(KIND=8) :: eqgrid(Ned,Npr,2) !grid for EoS REAL(KIND=8) :: x, y !grid point variables REAL(KIND=8) :: sp(2),ep(2) !start and end points REAL(KIND=8) :: pt1(2),pt2(2) !points for candidate generation REAL(KIND=8) :: expmin,expmax,expinc !exponent for grid generation REAL(KIND=8) :: rn !aux random number variable INTEGER(KIND=4) :: ydim(Ned),xdim !grid dimension INTEGER(KIND=4) :: candidate(NE,Ned) !EoS candidate LOGICAL :: ACCEPT !accept points for EoS LOGICAL :: GRIDMASK(Ned,Npr) !mask with upper and lower boundary encoded LOGICAL :: EOSMASK(Ned,Npr) !mask used throughout generation of EoS candidate INTEGER(KIND=4) :: seed !seed for random number generator CHARACTER(LEN=100) :: cdate !date used for seed CHARACTER(LEN=10) :: ctime !time used for seed INTEGER(KIND=4) :: i,j,l,m,n,lo,hi,idx,z !aux ints character (len=90) :: EOSfilename integer :: jfile,counter,file !INITIALIZATIONS !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !init random number generator with a date & time based seed cdate="" ctime="" CALL date_and_time(cdate,ctime) READ( ctime, '(I6)' ) seed CALL srand(seed) !Select upper band from: Sammarruca (NLO:case 1; N2LO: case 3; N3LO: case 5), Holt (NLO: case 7; N3LO: case 9; Jinniu (NLO; case 13; N3LO; case 15; N4LO: case 17) and Drischler (case 21) ! In all the files first column E in MeV/fm³; second column P in MeV/fm³ write(*,*) 'Select up band:1(SammaNLO);3(SammaN2LO);5(SammaN3LO);7(HoltNLO);9(HoltN3LO);13(HuNLO);15(HuN3LO);17(HuN4LO);21(Dris)' read (*,*) file call select_case (file) do i=1,Nd read (40,*) edensity(i), Bandroof(i) E(i) = edensity(i) enddo close(40) E(1) = edensity(1) Pmax= Bandroof(Nd) Emax= edensity(Nd) !Select lower band write(*,*) 'Select down band case: 2(SammaNLO);4(SammaN2LO);6(SammaN3LO);8(HoltNLO);& 10(HoltN3LO);14(HuNLO);16(HuN3LO);18(HuN4LO);22(Dris)' read (*,*) file call select_case (file) do i=1,Nd read (40,*) edensity (i), Bandground(i) enddo close (40) P(1) = Bandground(1) Pmin= Bandground(Nd) ! INITIALIZATION FINISHED, START COMPUTING PRESSURES !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccEOSSammarrucaNLOcutoff1GeV do jfile = 50+1,50+Nep !the 30 is to avoid low reserved numbers (such as the screen...) \ write (EOSfilename, '( "EOSDATA.DIR/EOSSammarrucaNLOcutoff1GeV", I6.6, ".dat" )' ) jfile-50 OPEN(unit=jfile,file=EOSfilename) enddo do m=1, Nep ! This loop generates each EoS one at a time 550 jfile=m+50 Rewind (jfile) !EOS points low band factor=Rand() do i=2,Nd bandfloor (i) = Min(bandground(i),P(i-1)) 150 P(i)= factor*bandroof(i)+(1-factor)*bandfloor(i) IF (P(i) .GT. bandroof(i) .OR. (P(i) .LT. bandground(i))) THEN factor=Rand() goto 150 ELSE IF(P(i) .GE. P(i-1)) THEN write(jfile,*) E(i), P(i) ELSE factor=Rand() goto 150 endif endif enddo close(40) close (44) !EOS points intermediate band !store boundary points for easy use, don't change P1x=Emax P2x=Emax P1y=Pmin P2y=Pmax P1intx=P1y-Q1y+Q1x !bndry pt P1inty=P1y !bndry pt P2intx=Q2y-P2y+P2x !bndry pt P2inty=Q2y !bndry pt P1(1) = P1x; P1(2) = P1y P2(1) = P2x; P2(2) = P2y Q1(1) = Q1x; Q1(2) = Q1y Q2(1) = Q2x; Q2(2) = Q2y P1int(1) = P1intx; P1int(2) = P1inty P2int(1) = P2intx; P2int(2) = P2inty !print upper and lower bound as function of energy density (x) !and write it to a file. OPEN (UNIT=10,FILE='boundary_pts.dat',ACTION="WRITE",STATUS="REPLACE") DO i=1, Nbd x = P2x + REAL((i-1),8)/(REAL(Nbd-1,8))*(Q2x - P1x) ! PRINT *,Real((i-1),8), x, up_bd(x,P2,Q2,P2int), lo_bd(x,P1,Q1,P1int) WRITE( 10,'(F18.5 X F18.5 X F18.5)') x, up_bd(x,P2,Q2,P2int), lo_bd(x,P1,Q1,P1int) END DO CLOSE(10) !Construct the grid. !The variable ydim is an array that stores the number of available points in y-direction !(pressure) for a given x-point (energy density). Available means that is lies within the !boundaries derived from causality and monotony. Initially, it is put to zero. !The data type is integer 4. ydim = 0 !The grid for the equation candidates that is to be constructed is also initialized to zero. !Its dimensions are: Ned x Npr x 2. For every grid index i and j, there is an energy-density !value ( indices (i,j,1) ) and a pressure value ( indices (i,j,2) ). !The data type is real 8. eqgrid = 0D0 !This logical mask keeps track of grid-points that are within the desired boundary. !It is used internally for the generation of the grid. It has dimensions ! Ned x Npr, and only those entries which satisfy the boundary constraints are !put to T later. GRIDMASK = .FALSE. !Next, fill the (logarithmic) grid with values that lie in the range !as given by the boundary points hardcoded above. DO i=1, Ned !loop over energy density points (x-axis) expmin = LOG10(P1x+50) expmax = LOG10(Q2x-1700) expinc = expmin + REAL((i-1),8)/(REAL(Ned-1,8))*(expmax-expmin) x = 10.0_8**(expinc) DO j=1,Npr !loop over pressure points (y-axis) expmin = LOG10(P1y+0.5) expmax = LOG10(Q2y-500) expinc = expmin + REAL((j-1),8)/(REAL(Ned-1,8))*(expmax-expmin) y = 10.0_8**(expinc) !PRINT *, x, y eqgrid(i,j,1) = x eqgrid(i,j,2) = y END DO END DO !Now fill the logical mask for all values that lie within the boundaries. DO i=1, Ned DO j=1, Npr x = eqgrid(i,j,1) y = eqgrid(i,j,2) IF( y .GT. lo_bd(x,P1,Q1,P1int) .AND. y .LT. up_bd(x,P2,Q2,P2int)) THEN GRIDMASK(i,j) = .TRUE. END IF END DO END DO !Write the grid data to file. OPEN (UNIT=20,FILE='EoS_grid.dat',ACTION="WRITE",STATUS="REPLACE") DO i=1, Ned DO j=1, Npr IF( GRIDMASK(i,j) .EQV. .TRUE. ) THEN WRITE (20,'(F18.5 X F18.5)') eqgrid(i,j,1), eqgrid(i,j,2) END IF END DO END DO CLOSE(20) !Determine number of x-points (it general, this should match Ned). xdim = 0 ydim = 0 ydim = COUNT(GRIDMASK,2) ! PRINT *, ydim DO i=1, Ned !Count all e-density values for which a non-zero number of pressure points satisfy the b.c. IF( ydim(i) .GT. 0 ) xdim = xdim + 1 END DO !! PRINT *, "Number of x-points: ", xdim !DO i=1, Ned ! PRINT *, "Number of y-points at xpos = ", i, " : ", ydim(i) !END DO !PRINT *, ydim !PRINT *, eqgrid(1,1,1), eqgrid(1,1,2) ! PRINT *, "Done." ! PRINT *, "=================================================" ! PRINT *, "Generating EoSs." !Now that the grid is set, let us create the EoS candidates sp(1)=P1x sp(2)=P(Nd) ! l=1 !This parameter l numbers the candidates from l=1,...,NE l=1 DO WHILE( l .LE. NE ) !loop over number of candidates EOSMASK = GRIDMASK !initialize aux logical mask ACCEPT = .FALSE. !upper boundary condition flag !come up with a random starting point (boundary condition) sp(1) = P1x sp(2) = P(Nd) !come up with a random end point (boundary condition) ep(1) = Q1x ep(2) =RAND()*(Q2y-Q1y)+Q1y !initialze aux ints i = 1 lo = 0 hi = 0 idx = 0 DO WHILE (i .LE. Ned+1) !loop over points in x-direction (energy density) !For a given position in energy density (x-direction), starting from !the randomly chosen starting point (boundary condition), determine !a possible paths forward by chosing points at random. If a chosen random !point (a random step) satisfies the monoyony and causality constraints, !accept the point and move on. If it does not, disregard the candidate and !start over. One could optimize this procedure by taking only one step back !and alter the last step choice, thereby systematically excluding candidates !that are to be excluded. However, even without this improvement the code !seems to run reasonably fast. ! print*, 'llegué al if ' IF( i .EQ. 1 ) THEN !i=1 corresponds to the (randomly chosen) start point pt1 = sp do j=1,ydim(i) pt2 = eqgrid(i,j,:) IF( check_constraints(pt1, pt2) .EQV. .FALSE. ) THEN EOSMASK(i,j) = .FALSE. ELSE IF( lo .EQ. 0 ) THEN lo = j END IF hi = j ENDIF ! PRINT *, "Accepted!" END DO ELSE IF( i .LE. Ned) THEN !i <=Ned are the points that are within the two boundary points. CALL RANDOM_NUMBER(rn) idx = FLOOR((hi+1-lo)*rn) + lo IF( MOD(l,2) .EQ. 0D0 ) THEN !This introduces the possibility to intorduce a bias for the random procedure. IF( idx .GT. lo+1 ) idx = FLOOR((hi+1-lo)*rn) + lo END IF lo = idx pt1 = eqgrid(i-1,idx,:) DO j=1, ydim(i) pt2 = eqgrid(i,j,:) IF( check_constraints(pt1, pt2) .EQV. .FALSE. ) THEN EOSMASK(i,j) = .FALSE. ELSE hi = j ENDIF END DO ELSE IF( i .EQ. Ned+1) THEN !Finally, check that also the last step (the choice of the end point) is OK. CALL RANDOM_NUMBER(rn) idx = FLOOR((hi+1-lo)*rn) + lo pt1 = eqgrid(i-1,idx,:) pt2 = ep IF( check_constraints(pt1, pt2) .EQV. .TRUE. ) THEN ACCEPT = .TRUE. ! PRINT *, "Accepted!" END IF END IF !Up to this point, we have determined wether there is a path forward for the chosen !random points. !PRINT *, "i = ", i !PRINT *, "lo = ", lo !PRINT *, "hi = ", hi !PRINT *, "idx = ", idx !PRINT *, "-------------------------------" !If the step above has been found to satisfy the constraints, accept the point (for now). IF ( ( (i .LE. Ned ) .AND. (COUNT(EOSMASK(i,:)) .GT. 0 ) ) .OR. (ACCEPT .EQV. .TRUE.) ) THEN ! print*, 'Paso por la primera rama ' IF ( i .GT. 1 ) THEN candidate(l,i-1) = idx END IF i = i+1 ACCEPT = .FALSE. ELSE !Else kicks in if the above procedure failed to produce a candidate. Reset and restart. ! print*, 'Paso por la segunda rama ' ! print*, 'comprobando bucle en i ',i, Ned+1 EOSMASK = GRIDMASK !initialize logical mask ACCEPT = .FALSE. !upper boundary condition flag ! REWIND(jfile) counter=counter+1 if (counter.LT.500) then sp(1) = P1x sp(2) = P(Nd) ep(1) = Q1x ep(2) = RAND()*(Q2y-Q1y)+Q1y i=1 lo = 0 hi = 0 idx = 0 else counter=0 Rewind (jfile) go to 550 end if END IF END DO !This point is only reached if a candidate has been produced successfully. !The coordinates of the candidate are written to a file. IF ( l .EQ. 1 ) THEN OPEN (UNIT=30,FILE="EoS_pts.dat",ACTION="WRITE",STATUS="REPLACE") WRITE(30,*) WRITE (30,'(F18.5 )',ADVANCE='NO') sp(1) DO i=1, xdim WRITE (30,'(F18.5 )',ADVANCE='NO') eqgrid(i,1,1) END DO WRITE (30,'(F18.5 )',ADVANCE='NO') ep(1) WRITE(30,*) WRITE (30,'(F18.5 )',ADVANCE='NO') sp(2) DO i=1, xdim WRITE (30,'(F18.5 )',ADVANCE='NO') eqgrid(i,candidate(l,i),2) write(jfile,*) eqgrid(i,1,1), eqgrid(i,candidate(l,i),2) END DO WRITE (30,'(F18.5 )',ADVANCE='NO') ep(2) CLOSE(30) ELSE OPEN (UNIT=30,FILE="EoS_pts.dat",ACTION="WRITE",POSITION="APPEND",STATUS="OLD") WRITE (30,'(F18.5 )',ADVANCE='NO') sp(2) DO i=1, xdim WRITE (30,'(F18.5 )',ADVANCE='NO') eqgrid(i,candidate(l,i),2) END DO WRITE (30,'(F18.5 )',ADVANCE='NO') ep(2) WRITE(30,*) CLOSE(30) END IF ! PRINT *, "l : ", l l = l + 1 !Increment the candidate id. ! print*, E(i),P(i) ! PRINT *, "Done." ! PRINT *, "=================================================" ! PRINT *, "Output files:" ! PRINT *, "=================================================" ! PRINT *, "boundary_pts.dat (contains upper and lower bdry)" ! PRINT *, "-------------------------------------------------" ! PRINT *, "EoS_grid.dat (contains grid points)" ! PRINT *, "-------------------------------------------------" ! PRINT *, "EoS_pts.dat (contains coordinates of candidates: " ! PRINT *, " 1st line: x-coords, rest: y-coords, " ! PRINT *, " each line corresponds to 1 candidate)" ! PRINT *, "=================================================" ! PRINT *, "Done." ! PRINT *, "=================================================" !EOS points high band open(unit=43,file='tabla_High.dat') !pQCD table at muB=2,6 GeV in three columns (E, Pmin, Pmax) in MeV/fm³ !open (unit=43,file='tabla_High1.dat') !pQCD table at muB=2.8 GeV in three columns (E, Pmin, Pmax) in MeV/fm³ do i=1,Neh read (43,*) table(i,1), table(i,2), table(i,3) ! if (stat /= 0) exit Ehigh(i) = table(i,1) Bandhighground(i) = table(i,2) Ehigh(1) = Q1x Bandhighroof(i) = table(i,3) Phigh(1) = ep(2) enddo write(jfile,*) Ehigh(1), Phigh(1) do i=2,Neh Bandhighfloor(i) = Min(bandhighground(i),Phigh(i-1)) 350 Phigh(i)= factor*bandhighroof(i)+(1-factor)*Bandhighfloor(i) IF(Phigh(i) .GT. bandhighroof(i) .OR.(Phigh(i) .LE. bandhighground(i))) THEN factor=Rand() goto 350 Else IF(Phigh(i) .GE. Phigh(i-1)) THEN write(jfile,*) Ehigh(i), Phigh(i) Else factor=Rand() goto 350 ! print*, rn endif endif enddo Close (43) enddo Close (jfile) Close (30) !Close (45) ! print*, 'aqui ', m enddo !finish the loop iterating over the equation of state END PROGRAM main