c*******Program 'den_147.f'
c*******Plot the surface-density contour map
c******* subroutine contour(a,idim,jdim,i1,i2,j1,j2,c,nc,tr)
c******* integer idim,jdim,i1,j1,i2,j2,nc
c******* real a(idim,jdim),c(*),tr(6)
c parameter(n1=147,nc=10)
parameter(n1=147,nc=6)
c parameter(n1=76,n3=5,n4=45,nc=10)
parameter(n_x=480,n_y=480)
c parameter(r0=3.0,n_win=15)
parameter(r0=4.0,n_win=20)
c*******sigma=0.5*n_win*********************
integer itype(n1),n_win
real ra(n1),dec(n1),x(n1),y(n1)
real xx(1),yy(1)
real x1(2),y1(2),x2(2),y2(2)
real tr(6),c(nc),a(n_x,n_y)
real aa(n_x,n_y)
c real xb3(n3),yb3(n3),czb3(n3)
c real xb4(n4),yb4(n4),czb4(n4)
c integer num3(n3),num4(n4)
c data
c/0.09,0.15,0.21,0.27,0.33,0.39,0.45,0.51,0.63,0.75/
c data c/0.09,0.15,0.21,0.27,0.33,0.39/
data c/0.05,0.09,0.13,0.17,0.21,0.25/
¡¡
c****************************************************************
c******* data tr/tr1,tr2,tr3,tr4,tr5,tr6/
c******* ===> data tr/-a,tr2,0.0,-a,0.0,tr6/
c******* ===> tr2=tr6=2a/n_x=2a/n_y
c***************************************************************
data tr/-60.0,0.25,0.0,-60.0,0.0,0.25/
open(1,file='x_y_cz_147.cat',
status='old')
c open(3,file='x_y_cz_fore.dat',status='old')
c open(4,file='x_y_cz_back.dat',status='old')
c open(2,file='zzzz')
c*******center of Abell 779 is ra<9:19:50.8> dec<33:46:17>
call pgbegin(0,'?',1,1)
call pgpaper(0.0,1.0)
call pgscf(2)
call pgslw(2)
call pgenv(80.0,-80.0,-80.0,80.0,0,0)
call pglabel('R.A. Offset (arcmin)','DEC.
Offset (arcmin)',
+
'Spatial distribution of A779 147 member galaxies')
call pgsci(1)
call pgslw(3)
call pgarro(65.0,65.0,65.0,75.0)
call pgarro(65.0,65.0,75.0,65.0)
call pgptxt(63.0,75.0,0.0,0.0,'N')
call pgptxt(76.0,58.0,0.0,0.0,'E')
call pgptxt(-50.0,74.0,0.0,0.0,'1 Mpc')
call pgmove(-40.0,72.0)
call pgdraw(-76.7,72.0)
call pgmove(-40.0,71.0)
call pgdraw(-40.0,73.0)
call pgmove(-76.7,71.0)
call pgdraw(-76.7,73.0)
call pgsci(9)
call pgsls(2)
call pgsfs(2) !good for plot circle
call pgcirc(0.0,0.0,78.0)
c call circ(0.0,0.0,78.0)
i=0
do 100 i=1,n1
read(1,*) nn,x(i),y(i),cz
100 continue
call pgpoint(n1,x,y,3)
call pgsls(1)
call
surf_density(n1,x,y,r0,n_x,n_y,a)
call smooth_gauss(n_win,n_x,n_y,a,aa)
z_max=0.0
z_min=100.0
do ii=1,n_x
do jj=1,n_y
c write(2,'(f10.6)')
aa(ii,jj)
if(aa(ii,jj).gt.z_max)
z_max=aa(ii,jj)
if(aa(ii,jj).lt.z_min)
z_min=aa(ii,jj)
end do
end do
write(*,*) z_min,z_max
xx(1)=0.0
yy(1)=0.0
call pgslw(4)
call pgsch(3.0)
call pgpt(1,xx,yy,2)
call pgslw(2)
call
pgcont(aa,n_x,n_y,1,n_x,1,n_y,c,nc,tr)
call pgend
888 format(i4,2(1x,f11.6),1x,f8.1)
close(1)
c close(2)
c close(3)
c close(4)
end
subroutine smooth(n_win,n_x,n_y,z,zz)
integer n_win,n_x,n_y
real z(n_x,n_y),zz(n_x,n_y)
do i=1,n_x
do j=1,n_y
zz(i,j)=0.0
end do
end do
do 100 i=1,n_x
do 200
j=1,n_y
k1=kmax(1,i-(n_win-1)/2)
k2=kmin(n_x,i+(n_win-1)/2)
m1=kmax(1,j-(n_win-1)/2)
m2=kmin(n_y,j+(n_win-1)/2)
i_count=0
do 300 k=k1,k2
do 400
m=m1,m2
dist=sqrt((k-i)*(k-i)*1.0+(m-j)*(m-j)*1.0)
if(dist.le.n_win) then
i_count=i_count+1
zz(i,j)=zz(i,j)+z(k,m)
end if
400 continue
300 continue
zz(i,j)=zz(i,j)/i_count
200 continue
100 continue
return
end
function kmax(n1,n2)
integer n1,n2
kmax=n2
if(n1.ge.n2) kmax=n1
return
end
function kmin(n1,n2)
integer n1,n2
kmin=n2
if(n1.le.n2) kmin=n1
return
end
subroutine
surf_density(n1,x,y,r0,n_x,n_y,z)
integer n1,n_x,n_y
real x(n1),y(n1),r0
real z(n_x,n_y)
real dist,f
f=n_x/120.0
pi=3.1415926535
area=pi*r0*r0
do i=1,n_x
do j=1,n_y
z(i,j)=0.0
end do
end do
do i=1,n1
x(i)=x(i)+60.0
y(i)=y(i)+60.0
end do
do 100 i=1,n_x
do 200
j=1,n_y
do 300 k=1,n1
dist=sqrt((x(k)*f-i)**2.0+(y(k)*f-j)**2.0)
c
dist=sqrt((x(k)*2.0-i/10.0)**2.0+(y(k)*2.0-j/10.0)**2.0)
if(dist.gt.r0*f) goto 300
z(i,j)=z(i,j)+1.0/area
300 continue
200 continue
100 continue
return
end
function gauss1(x,sigma)
real x, sigma,pi,pp
pi=3.1415926535
pp=exp(-1.0*x*x/(2.0*sigma*sigma))
gauss1=pp/sqrt(2.*pi)/sigma
return
end
subroutine
smooth_gauss(n_win,n_x,n_y,z,zz)
integer n_win,n_x,n_y
real z(n_x,n_y),zz(n_x,n_y)
real frac,all_frac,sigma
sigma=0.5*n_win
do i=1,n_x
do j=1,n_y
zz(i,j)=0.0
end do
end do
do 100 i=1,n_x
do 200
j=1,n_y
k1=kmax(1,i-(n_win-1))
k2=kmin(n_x,i+(n_win-1))
m1=kmax(1,j-(n_win-1))
m2=kmin(n_y,j+(n_win-1))
i_count=0
all_frac=0.0
do 300 k=k1,k2
do 400 m=m1,m2
dist=sqrt((k-i)*(k-i)*1.0+(m-j)*(m-j)*1.0)
all_frac=all_frac+gauss1(dist,sigma)
400 continue
300 continue
do 301 k=k1,k2
do 401
m=m1,m2
dist=sqrt((k-i)*(k-i)*1.0+(m-j)*(m-j)*1.0)
frac=gauss1(dist,sigma)/all_frac
zz(i,j)=zz(i,j)+frac*z(k,m)
401 continue
301 continue
c
zz(i,j)=zz(i,j)/i_count
200 continue
100 continue
return
end
subroutine circ(x,y,r)
real x,y,r
real xx(500),yy(500)
pi=3.1415926535
d_theta=2*pi/500.
do i=1,500
xx(i)=x+r*cos(i*d_theta)
yy(i)=y+r*sin(i*d_theta)
end do
call pgline(500,xx,yy)
return
end
back |