# Script to produce simple plots to demonstrate how stacking spectra improves sensitivity 
# Produces an output file containing all individual spectra plus their stack

import numpy

speclength = 100	# Number of elements in each spectrum
nspec  = 100		# Number of spectra to generate and plot
sics   = 45			# Start channel to inject signal
sicf   = 55			# Final channel to inject signal
sigmai = 1.0		# Signal to inject between sics and sicf

# Array containing all the spectra. The stacked spectrum will be computed afterwards
mainspec = numpy.zeros((speclength, nspec))	# rows, columns

for i in range(nspec):
	# Generate spectra of Gaussian noise. Loc = mean, scale = standard deviation, size = number of elements
	newspec = numpy.random.normal(loc=0.0, scale=1.0, size=speclength)
	mainspec[:,i] = newspec[:]
	
	# Inject a signal into the specified channel range
	mainspec[45:55,i] = mainspec[sics:sicf,i] + sigmai
	
	
# Do the stacking across each ROW, saving the final spectrum as its own, new array
stackplot = []
for i in range(0, speclength):
	stacked = numpy.sum(mainspec[i,:]) / float(nspec)
	stackplot.append([i, stacked])
	
stackplot = numpy.asarray(stackplot)


# Save the indicies, all individual spectra, and stacked spectrum to a file	
stackedplotfile = open('StackingPlotAllSpec_'+str(speclength)+'.txt','w')

# Give it a header
headerstart = '#Index,'
headerend   = 'Stacked'
headermid   = ''

for i in range(nspec):
	headermid = headermid+'Spectra_'+str(i)+','

stackedplotfile.write(headerstart+headermid+headerend+'\n')

for i in range(len(stackplot)):
	individualspec = ''
	for j in range(speclength):
		individualspec = individualspec+str(mainspec[i,j])+','
		
	stackedplotfile.write(str(i)+','+individualspec+str(stackplot[i,1])+'\n')
	
stackedplotfile.close()