BubbleExoplanets.r

BubbleExoplanets.r

# source("C:/Users/LucaF/Documents/My Dropbox/Software/R/Exoplanets.r")
rm(list=ls(all=TRUE))

if (!exists('palette_l')) {palette_l<-50} # how many steps in the palette

N<-650 # number of exoplanets
R<-1000 # radius of a 2D-circle where you'll want to plot them
MinMass<-5; # minimum mass for the planets
MaxMass<-100; # maximum mass for the planets
P<-6 # power exponent for the mass distribution
Q<-3 # this exponent scales the mass to radii...
margin<-mean(c(MinMass,MaxMass)); # a minimum separation you want to mantain between planets when plotted
# a function to name each planet with a random string of letters (optional)
getRandString<-function(len=8) return(paste(sample(c(LETTERS,letters),len,replace=TRUE),collapse=''))

exoplanets<-list(Mass=sort(runif(N, MinMass^(1/P), MaxMass^(1/P))^P, decreasing=TRUE), Name=replicate(n=N, getRandString()))

mtrx<-matrix(nrow=N, ncol=2); colnames(mtrx)<-c('x','y')

getxy_withinR<-function(R=10) {
	r<-runif(1,0,R^2)^(1/2);	theta<-runif(1,0,2*pi);
	x<-r*cos(theta)
	y<-r*sin(theta)
	return(c(x,y))
}

# stop("It's OK")

check_bounce <- function(x,y,rad,X,Y,RAD) {
		distance<-sqrt((X-x)^2 + (Y-y)^2); # print (paste("distance is ",distance))
		radsum<-RAD^(1/Q)+rad^(1/Q)+margin; # print(paste("radius sum is ",radsum))
		return(sum(distance < radsum))
}

lim<-R+MaxMass
plot(NULL, xlim=c(-lim,lim), ylim=c(-lim,lim))

for (n in 1:N) {
	print(paste("placing",n, exoplanets$Name[n]))
	# generate new coordinates 
	xy<-getxy_withinR(R); x<-xy[1]; y<-xy[2];
	if (n!=1) {
		tocheck<-(1:(n-1))
		while(overlap<-check_bounce(x,y,exoplanets$Mass[n],mtrx[tocheck,'x'],mtrx[tocheck,'y'],exoplanets$Mass[tocheck]) > 0)
		{
			xy<-getxy_withinR(R); x<-xy[1]; y<-xy[2];
		}
	}
	mtrx[n,'x']<-x;	mtrx[n,'y']<-y;
	# draw a circle once you know where
	exoplanets$X[n]<-xy[1]; exoplanets$Y[n]<-xy[2]; 
	draw.circle(x=exoplanets$X[n], y=exoplanets$Y[n], r=exoplanets$Mass[n]^(2/Q), col=rgb(runif(1,0,1),runif(1,0,1),runif(1,0,1),(0.5+n/(2*N))), border=rgb(runif(1,0,1),runif(1,0,1),runif(1,0,1),(0.5+n/(2*N))))
	# Sys.sleep(0.03)
}

exoplot<-ggplot(as.data.frame(exoplanets))+ geom_point(aes(x=X, y=Y, size=Mass, color=Name))+ scale_size(range = c(MinMass, MaxMass)) +opts(legend.position='none')#+coord_polar()
# print(exoplot) # doesn't look the way I want it. Pity.