public
Last active

Observing Dark Worlds visualization

  • Download Gist
plot_skies.R
R
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
################## Plot training skies ###################
##
## corey.chivers@mail.mcgill.ca
##
##########################################################
 
## calculate a vector given
## x,y,e1,e2
gal_line<-function(g,scale=100)
{
# e1 is stretch along 0degrees (+1)
# or along 90degrees (-1)
# e2 is stretch along 45 (+1)
# or along 135 (-1)
 
#e1
if(g$e1<=0)
{
x<-0
y<-scale*g$e1
}
if(g$e1>=0)
{
x<-scale*g$e1
y<-0
}
 
#e2
x<-x+scale*g$e2*cos(pi/4)
y<-y+scale*g$e2*sin(pi/4)
#start and end points
l_seg<-rbind(c(g$x+x,g$y+y),c(g$x-x,g$y-y))
 
return(l_seg)
}
add_halo<-function(dm,n_rings=20,max_ringsize=2000)
{
d<-seq(0.1,1,length.out=n_rings)
radii<-1/d
radii<-radii/max(radii) * max_ringsize
for(i in 1:dm$numberHalos)
{
centre<-c(dm[1,2*(i-1)+5],dm[1,2*(i-1)+6])
for(r in radii)
{
x=r*cos(seq(0,2*pi,length.out=100)) + centre[1]
y=r*sin(seq(0,2*pi,length.out=100)) + centre[2]
polygon(x,y,col=rgb(0,0,1,0.05),border=NA)
}
}
}
 
plt_sky<-function(sky,dm,...)
{
## Plot the galaxies
par(bg='black',col='white',col.main='white')
 
plot(sky$x,sky$y,col='black',...)
 
## Add the Dark Matter
add_halo(dm)
 
for(i in 1:nrow(sky))
{
gl<-gal_line(sky[i,])
lines(gl,col='white')
#print(dist(gl))
}
points(sky$x,sky$y,col='grey',pch=20,cex=0.5)
}
##########################################################################
 
 
tr_halos<-read.csv('Training_halos.csv')
tr_folder<-'Train_Skies/'
tr_files<-paste(tr_folder,list.files(tr_folder),sep='')
 
## Extract sky number from file naming convention
index_sky<-1:length(tr_files)
for(i in 1:length(tr_files))
{
m <- regexpr('[0-9]+', tr_files[i])
index_sky[i]<-as.integer(regmatches(tr_files[i], m))
}
tr_files<-tr_files[order(index_sky)]
####
 
 
pdf('plots/training_skies.pdf')
for(i in 1:length(tr_files))
{
sky<-read.csv(tr_files[i])
dm<-tr_halos[i,]
plt_sky(sky,dm,main=paste('Sky',i))
print(paste('plotting sky',i))
}
dev.off()

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.