User:Zeracles/Background/Code/Stars

See also this quick download

%background.m
%by Anthony Smith (Zeracles); astroant@hotmail.com, asmith@physics.usyd.edu.au, http://www.physics.usyd.edu.au/~asmith/, http://wiki.uqm.stack.nl/User:Zeracles
%see http://www.star-control.com/forum/index.php/topic,1660.0.html
%nice stellar geometry can be obtained with anglemagnificationbase=0.5;anglemagnificationindex=5;anglecutoffradiusfactor=1;

n=800;
resolution=1;

xmin=0;
xmax=1600;
ymin=0;
ymax=1200;

xextent=xmax-xmin;
yextent=ymax-ymin;

numxpix=round((xextent)/resolution);
numypix=round((yextent)/resolution);

spokethreshold=0;
spokes=4;
anglemagnificationbase=0.5;
anglemagnificationindex=5;

cutoffradiusfactor=4;
anglecutoffradiusfactor=1;

plotposns=1;
plotimage=1;

xposns=xmin+(xmax-xmin)*rand(n,1);
yposns=ymin+(ymax-ymin)*rand(n,1);

if plotposns==1
	figure
	box on
	plot(xposns,yposns,'k.')
end

redmap=zeros(numypix,numxpix);
greenmap=zeros(numypix,numxpix);
bluemap=zeros(numypix,numxpix);

for i=linspace(1,n,n)
	  i

	  xposn=xposns(i);
	  yposn=yposns(i);

	  radius=((exp(rand(1)))*14-0.01)/2;
	  intensity=2*rand(1);

	  redness=rand(1);
	  greenness=rand(1);
	  blueness=rand(1);
	
	  spokeprobability=rand(1);

	  spokeswitch=spokeprobability>spokethreshold;

	  spokeangleoffset=pi/4;
	
        ygaussiancontribution=[];
        for gridxpixel=linspace(1,numxpix,numxpix)
                gridxposn=xmin+(gridxpixel-1)*(xextent)/(numxpix-1);
                xgaussiancontribution=[];
                for gridypixel=linspace(1,numypix,numypix)
                      gridyposn=ymin+(gridypixel-1)*(yextent)/(numypix-1);
				
			xdiff=xposn-gridxposn;
			ydiff=yposn-gridyposn;
				
                      distance=((xdiff)^2+(ydiff)^2)^(1/2);
			distancefactor=distance/radius;
			if distancefactor>cutoffradiusfactor
				littlexgaussiancontribution=0;
			end
			if distancefactor<=cutoffradiusfactor
                                littlexgaussiancontribution=exp(-1*(distance^2)/(2*radius^2));
					if spokeswitch==1
						anglesweep=2*pi/spokes;
						angle=atan(ydiff/xdiff);
						correctedangle=angle-spokeangleoffset;
						angledisplacementnumber=round((correctedangle)/(anglesweep));
						angledisplacement=abs((correctedangle)-angledisplacementnumber*(anglesweep));
						anglemagnification=anglemagnificationbase+(1-anglemagnificationbase)*((1-anglecutoffradiusfactor*distancefactor*(angledisplacement*(anglesweep/2))))^anglemagnificationindex;
						if anglemagnification<0
							anglemagnification=0;
						end
						littlexgaussiancontribution=littlexgaussiancontribution*anglemagnification;
					end
			end
                      xgaussiancontribution=[xgaussiancontribution;littlexgaussiancontribution];
                end
                ygaussiancontribution=[ygaussiancontribution,xgaussiancontribution];
        end
        gaussiancontribution=intensity*ygaussiancontribution;

	  redmap=redmap+gaussiancontribution*redness;
        greenmap=greenmap+gaussiancontribution*greenness;
        bluemap=bluemap+gaussiancontribution*blueness;
end

excessmask=redmap>1;
redmap=redmap-excessmask.*(redmap-1);

excessmask=greenmap>1;
greenmap=greenmap-excessmask.*(greenmap-1);

excessmask=bluemap>1;
bluemap=bluemap-excessmask.*(bluemap-1);

backgroundmap=cat(3,redmap,greenmap,bluemap);

if plotimage==1
	figure
	box on
	image(backgroundmap)
end