Code:
program grprofilebasic
implicit none
character*200 path,infile,pathh,hfile
character*8 sim, icfile
character*3 snap
integer*4 fn,nin,iargc, i,j
integer*4 nfile, nh,ng,ng200,nh200, nc,nall1
double precision rsep, dlgr,zred,ared,omegamz,rhobgdm,rhocritdm,rhocritg, rhobg,rhocritz,fgas
double precision, parameter :: omegam0 = 2.5d-1
double precision, parameter :: omegab = 4.5d-2
double precision, parameter :: h = 7.3d-1
double precision, parameter :: lgrmax = 0.0d0
double precision, parameter :: lgrmin = -1.3d0 ! r/r200 = 0.05
double precision, parameter :: rmax = 1.d0
double precision, parameter :: rmin = 5.d-2
integer, parameter :: nbin = 19
double precision lgr(nbin),tbin(nbin),lgrhodm(nbin),lgrhog(nbin),volmpc(nbin),&
volbin(nbin),hy(nbin),dytemp(nbin),ytemp(nbin),rtemp,r(nbin),drbin,
xfit(nbin)
integer etag(nbin),etadm(nbin), ibin
double precision yint,qsum,lgrhotot(nbin),sig(nbin),rmean(nbin),lgrmean(nbin)
integer*2, allocatable :: veldmin(:,:),velgasin(:,:),gasin(:,:)
real*4, allocatable :: posdm(:,:),posgas(:,:),temp(:),u(:),rg(:),rh(:),rho(:),hs(:), veldm(:,:)
!--- recall gasin(1) is u; gasin(2) is rho, gasin(3) is hsml
integer*4, allocatable :: idm(:),idg(:)
real*4 cen(1:3)
double precision dx,dy, dz, delta_c, nfwdeltac
integer idelta_c
character*3 delta_cfile
integer*4 hrank
integer*4 satflag
real*4 m200,mg,r200
real, parameter :: rcore = .01 ! rc/r200
double precision thalo, elhalo, v2sum,v2bar,dvphys,com(1:3),sigmadm, vbar(1:3), tcore
!---------------------------------------
! Constants for conversion to cgs units
double precision, parameter :: mgaspart = (3.12d9/h) ! M_sun/h, not the internal m units
double precision, parameter :: mdmpart = (1.422d10/h)
double precision, parameter :: msung = 1.989d33
double precision, parameter :: mgasg = (mgaspart*msung)
double precision, parameter :: mp = 1.6726d-24 ! g
double precision, parameter :: mu = 5.9d-1
double precision, parameter :: unitmassg = 1.989d43 ! g/h, recall 10^10 M_sun/h
double precision, parameter :: unitelcm = 3.085678d24 !
double precision, parameter :: unitvelcms = 1.d5
double precision, parameter :: unittimes = (unitelcm/unitvelcms)
double precision, parameter :: unitecgs = (unitmassg*(unitelcm**2)/(unittimes)**2)
double precision, parameter :: gamma = (5.d0/3.d0)
double precision, parameter :: kb = (1.3806d-16) ! boltzmann const in cgs
double precision, parameter :: unittempk = (((gamma-1.d0)*mu*mp/kb)*unitecgs/unitmassg)
double precision, parameter :: kelvinkev = 8.62e-8 ! multiply to go from K to keV
double precision, parameter :: lambda0 = 1.2d-24 ! cgs units for Lamba(T)
double precision, parameter :: pi = 3.14159d0
!---------Conversion functions
real*4 convertu,convertrho,converths
dlgr = (lgrmax-lgrmin)/(dble(nbin))
do i = 1, nbin
lgr(i) = lgrmin + (dble(i)-1.d0)*dlgr
r(i) = rmin + (dble(i)-1.d0)*drbin
enddo
do i = 1, nbin
lgrhodm(i) = 0.
lgrhog(i) = 0.
lgrhotot(i) = 0.d0
sig(i) = 0.
if (i.lt.nbin) then
volbin(i) = (4.d0*pi/3.d0)*((10.d0**(lgr(i+1)))**3 - (10.d0**(lgr(i)))**3)
endif
enddo
volbin(nbin) = (4.d0*pi/3.d0)*((10.d0**(lgr(nbin)+dlgr))**3 - (10.d0**(lgr(nbin)))**3)
call getarg(1,sim)
call getarg(2,snap)
fn = 1 ! hrank
delta_c = 2.d2
idelta_c = int(delta_c)
write (delta_cfile,'(i3)') idelta_c
path = '/phys/ms1/' // sim(1:len_trim(sim)) // '/' // snap // '/grp/'
pathh = '/phys/ms1/' // sim(1:len_trim(sim)) // '/' // snap // '/'
hfile = pathh(1:len_trim(pathh)) // 'mill_' // sim(1:len_trim(sim)) // '_' // snap // '.header'
open (unit = 10,file = hfile)
read (10,*) nall1
read (10,*) ared,zred
omegamz = 1./(1.+((1./omegam0)-1.)*ared**3)
rhobgdm = 5.e8*mdmpart/((500./h)**3) ! (Msun/h) / (Mpc/h)**3
rhocritdm = rhobgdm/omegamz
rhocritg = rhocritdm*omegab
rhobg = (5.e8*mdmpart + 5.e8*mgaspart)/((500./h)**3)
rhocritz = (rhobg/omegamz)
write (icfile, '(i8)') fn
infile = path(1:len_trim(path)) // 'grpbulk_' // delta_cfile // '_' // snap // '.' // icfile(verify(icfile,' '):8)
open (unit = 20, file = infile, form='unformatted', status = 'old', err = 100)
do i = 1, nbin
tbin(i) = 0.
etadm(i) = 0
etag(i) = 0
lgrhodm(i) = 0.
lgrhotot(i) = 0.
lgrhog(i) = 0.
sig(i) = 0
rmean(i) = 0.d0
lgrmean(i) = 0.d0
enddo
nfile = fn
hrank = fn
read (20) nh
allocate(posdm(1:3,1:nh))
allocate(veldmin(1:3,1:nh))
allocate(veldm(1:3,1:nh))
allocate(idm(1:nh))
allocate(rh(1:nh))
read (20) m200
read (20) mg
read (20) r200
read (20) cen
read (20) posdm
read (20) veldmin
read (20) idm
!cengrp(1:3,fn) = cen(1:3)
!rgrp(fn) = r200
fgas = mg/m200
satflag = 0
do j = 1, nh
dx = abs(posdm(1,j)-cen(1)) ! deal w/periodic boundaries
if (dx.gt.(.5)) then
dx = 1.-dx
endif
dy = abs(posdm(2,j)-cen(2))
if (dy.gt.(.5)) then
dy = 1.-dy
endif
dz = abs(posdm(3,j)-cen(3))
if (dz.gt.(.5)) then
dz = 1.-dz
endif
rh(j) = sqrt(dx**2 + dy**2 + dz**2)
veldm(1:3,j) = float(veldmin(1:3,j))
enddo
i = 1
nh200 = 0
do while (rh(i).le.r200)
ibin = int((log10(rh(i)/r200)-lgrmin)/dlgr) + 1 ! determine which bin to put in
if ((ibin.le.nbin).and.(ibin.gt.0)) then
etadm(ibin) = etadm(ibin) + 1 ! counting particles in bin
lgrhodm(ibin) = lgrhodm(ibin) + 1. ! dark matter density
lgrhotot(ibin) = lgrhotot(ibin) + mdmpart ! total density
sig(ibin) = sig(ibin) + 1
nh200 = nh200 + 1
rmean(ibin) = rmean(ibin) + (rh(i)/r200)
endif
i = i + 1
enddo
read (20) ng
allocate(posgas(1:3,1:ng))
allocate(velgasin(1:3,1:ng))
allocate(idg(1:ng))
allocate(gasin(1:3,1:ng))
allocate(rg(1:ng))
read (20) posgas
do j = 1, ng
dx = abs(posgas(1,j)-cen(1)) ! deal w/periodic boundaries
if (dx.gt.(.5)) then
dx = 1.-dx
endif
dy = abs(posgas(2,j)-cen(2))
if (dy.gt.(.5)) then
dy = 1.-dy
endif
dz = abs(posgas(3,j)-cen(3))
if (dz.gt.(.5)) then
dz = 1.-dz
endif
rg(j) = sqrt(dx**2 + dy**2 + dz**2)
enddo
read (20) velgasin
read (20) idg
read (20) gasin
close (20)
allocate(temp(1:ng))
allocate(u(1:ng))
allocate(rho(1:ng))
allocate(hs(1:ng))
thalo = 0.d0
tcore = 0.d0
elhalo = 0.d0
i = 1
ng200 = 0
nc = 0
do while (rg(i).le.r200)
ibin = int((log10(rg(i)/r200)-lgrmin)/dlgr) + 1
if ((ibin.le.nbin).and.(ibin.gt.0)) then
etag(ibin) = etag(ibin) + 1
lgrhog(ibin) = lgrhog(ibin) + 1.
lgrhotot(ibin) = lgrhotot(ibin) + mgaspart
rmean(ibin) = rmean(ibin) + (rg(i)/r200)
sig(ibin) = sig(ibin) + 1
u(i) = convertu(gasin(1,i)) ! should give u in cgs
rho(i) = convertrho(gasin(2,i))
hs(i) = converths(gasin(3,i))
temp(i) = u(i)*unittempk*kelvinkev ! should give T in keV
thalo = thalo + dble(temp(i))
tbin(ibin) = tbin(ibin) + dble(temp(i))
elhalo = elhalo + (rho(i)*sqrt(temp(i)))*unitmassg/(unitelcm**3) ! cgs
ng200 = ng200 + 1
endif
i = i+1
enddo
elhalo = elhalo*lambda0*mgasg/(mu*mp)**2
do i = 1, nbin
if (etag(i).gt.0) then
tbin(i) = tbin(i)/dble(etag(i))
rmean(i) = rmean(i)/sig(i)
lgrmean(i) = log10(rmean(i))
volmpc(i) = ((5.d2*r200/h)**3)*volbin(i)
lgrhotot(i) = lgrhotot(i)/volmpc(i)/rhocritz ! (Mass/Vol)/rhocrit
sig(i) = 1./sqrt(float(sig(i)))
else
lgrmean(i) = lgr(i)
lgrhotot(i) = 0.
sig(i) = 1.
endif
enddo
deallocate(posdm)
deallocate(veldmin)
deallocate(idm)
deallocate(veldm)
deallocate(posgas)
deallocate(velgasin)
deallocate(idg)
deallocate(gasin)
deallocate(temp)
deallocate(u)
deallocate(rho)
deallocate(hs)
deallocate(rg)
deallocate(rh)
end program grprofilebasic
real function convertu(gasin)
!----------------------
! Converts int*2 u to real*4
!----------------------
implicit none
integer*2 gasin
real*4, parameter :: lnumax = 22.
real*4, parameter :: lnumin = 3.
real*4, parameter :: didlnu = (float(65535)/(lnumax-lnumin))
convertu = exp(lnumin + (float(gasin)+float(32768))/didlnu)
return
end function convertu
real function convertrho(gasin)
!----------------------------
! Converts int*2 rho to real*4
!-------------------------
implicit none
integer*2 gasin
real*4, parameter :: lnrhomax = 10.
real*4, parameter :: lnrhomin = -4.
real*4, parameter :: didlnrho = (float(65535)/(lnrhomax-lnrhomin))
convertrho = exp(lnrhomin + (float(gasin)+float(32768))/didlnrho)
return
end function convertrho
real function converths(gasin)
!-----------------------------
! Converts in*2 hsml to real*4
!------------------------
implicit none
integer*2 gasin
real*4, parameter :: hscale = 5.e3
converths = float(gasin)/hscale
return
end function converths
double precision function nfwdeltac(c)
!------------------------
! Calculates deltac(c) for NFW profile, for a Delta=200 halo
!------------------------------------------
double precision c
nfwdeltac = (2.d2/3.d0)*(c**3/(log(1.d0+c) - (c/(1.d0+c))))
return
end function nfwdeltac
subroutine funcs(x,a,y,dyda,ma)
integer ma
double precision x,a(ma),y,dyda(ma),y2,a2
double precision nfwdeltac
double precision, parameter :: eps = 1.d-2
y = nfwdeltac(a(1))/(a(1)*1.d1**(x))/(1.+a(1)*1.d1**(x))**2
!y = nfwdeltac(a(1))/(a(1)*x)/(1.+a(1)*x)**2
a2 = a(1) + eps
y2 = nfwdeltac(a2)/(a2*1.d1**(x))/(1.+a2*1.d1**(x))**2
!y2 = nfwdeltac(a2)/(a2*x)/(1.+a2*x)**2
dyda(1) = (y2-y)/(a2-a(1))
return
end subroutine funcs
------------------------------------------
END OF CODE
------------------------------------------PS.
I would like to understand more indepth, what is really going on; but that is not my main problem. I want to write to a file not within the path described by the 'path', but on my desktop. How do I do that? (I have less than the basic understanding of C, Java)
