import string, math
from math import *

def closedpoints(func, a, b, Wavelength, Mirror_diameter):         # f(x)=func(x)
    TOL=1e-4
    h = b - a
    old2 = old = h * (func(a, Wavelength, Mirror_diameter) + func(b, Wavelength, Mirror_diameter)) / 2.0
    count = 0
    while 1:
        h = h / 2.0
        x, sum = a + h, 0
        while x < b:
            sum = sum + func(x, Wavelength, Mirror_diameter)
            x = x + 2 * h
        new = old / 2.0 + h * sum
        new2 = (4 * new - old) / 3.0
        if abs(new2 - old2) < TOL * (1 + abs(old2)):
		return (new2)
        old = new       # Trapezoid
        old2 = new2     # Simpson
        count = count + 1
	#Temp='closedpoints(%d): Trapezoid=%s, Simpson=%s' % (count, new, new2)
	#print Temp

def bessel_low(x):
#Returns the bessel function J1(x) for any real x.
	ax=0.0
	z=0.0
	xx=0.0
	y=0.0
	ans=0.0
	ans1=0.0
	ans2=0.0 
	if ((ax==fabs(x)) < 8.0): #Direct rational approximation.
		y=x*x
		ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))))
		ans2=144725228442.0+y*(2300535178.0+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0))))
		ans=ans1/ans2
	else:
		z=8.0/ax
		y=z*z
		xx=ax-2.356194491
		ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4+y*(0.2457520174e-5+y*(-0.240337019e-6))))
		ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6)))
		ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2)
		if (x < 0.0):
			ans = -ans
	return(ans)
		
def bessel_high(x):
	#global e
	e = 0.0
	#global e3
	e3 = 1e-8
	#global J0
	J0=0.0
	#global J1
	J1=0.0
	#global n
	n = 0
	#global m1
	#global m2
	#global n1
	#global n2
	a = 1
	a1 = 1
	a2 = 1
	b = 1
	c = 1
	e1 = 1e6
	e2=0.0
	e4=0.0
	m = -1
	x1 = 1.0 / (8.0 * x)
	x1 = x1 * x1
	m1 = 1.0
	m2 = 1.0
	n1 = -1.0 / (8.0 * x)
	n2 = -3.0 * n1
#	print e1, e3
	while 1:#(e1 < e3):
		m = m + 2
		a = a * m * m
		m = m + 2
		a = a * m * m
		c = c * x1
		a1 = a1 * a2
		a2 = a2 + 1.0
		a1 = a1 * a2
		a2 = a2 + 1.0
		e2 = a * c / a1
		e4 = 1.0 + (m + 2) / m + (m + 2) * (m + 2) / (a2 * 8 * x) + (m + 2) * (m + 4) / (a2 * 8 * x)
		e4 = e4 * e2
		if (fabs(e4) > e1):
			break
		e1 = fabs(e2)
		m1 = m1 - e2
		m2 = m2 + e2 * (m + 2) / m
		n1 = n1 + e2 * (m + 2) * (m + 2) / (a2 * 8 * x)
		n2 = n2 - e2 * (m + 2) * (m + 4) / (a2 * 8 * x)
		n=n+1
		if (e1 < e3):
			break
	a=pi
	e = e2
	b = sqrt(2.0 / (a * x))
	#J0 = b * (m1 * cos(x - a / 4) - n1 * sin(x - a / 4.0))
	J1 = b * (m2 * cos(x - 3 * a / 4.0) - n2 * sin(x - 3 * a / 4.0))
	return(J1)

def bessel(r):
	if fabs(r)<1.5 and r != 0:
		Bessel=sign(r)*bessel_low(sign(r)*r)
	if fabs(r)>=1.5:
		Bessel=sign(r)*bessel_high(sign(r)*r)
	if r==0:
		Bessel=0.0
	return(Bessel)

def sign(x):
	if x<0:
		Temp=-1.0
	if x>=0:
		Temp=1.0
	return(Temp)

def air_func(x, y, Wavelength, Mirror_diameter):
	Parameter=(pi*Mirror_diameter/Wavelength) # Instrument parameter
	Arcsec_multiplier=(2.*pi/(3600.*360.))
	r=Parameter*Arcsec_multiplier*sqrt(x*x+y*y)
	Normal=pi*Mirror_diameter*Mirror_diameter/(4.*Wavelength*Wavelength)*Arcsec_multiplier*Arcsec_multiplier
	if x==0. and y==0.:
		Temp=1.
	else:
		Temp=(2.*bessel(r)/r)*(2.*bessel(r)/r)
	Temp=Temp*Normal
	return (Temp)

def air_func_int(Limit_x1, Limit_x2, Limit_y1, Limit_y2, Wavelength, Mirror_diameter):
	Step=abs(Limit_x1-Limit_x2)/100.
	x=Limit_x1
	y=Limit_y1
	Integral=0.
	while x<Limit_x2:
		while y<Limit_y2:
			Value=air_func(x, y, Wavelength, Mirror_diameter)
			Integral=Integral+Value*Step*Step
			#Temp="%s\t%s\t%s\n" % (x, y, Funktion)
			#File.write(Temp)
			y=y+Step
		y=Limit_y1
		x=x+Step
	return(Integral)


def air(r, Wavelength, Mirror_diameter):
	Parameter=(pi*Mirror_diameter/Wavelength) # Instrument parameter
	Arcsec_multiplier=(2.*pi/(3600.*360.))
	r=Parameter*Arcsec_multiplier*r
	Normal=pi*Mirror_diameter*Mirror_diameter/(4.*Wavelength*Wavelength)*Arcsec_multiplier*Arcsec_multiplier
	if r==0.:
		Temp=1.
	else:
		Temp=(2.*bessel(r)/r)*(2.*bessel(r)/r)
	Temp=Temp*Normal
	return (Temp)

def air_int(r, Wavelength, Mirror_diameter):
	Step=r/500.
	Parameter=(pi*Mirror_diameter/Wavelength) # Instrument parameter
	Arcsec_multiplier=(2.*pi/(3600.*360.))
	r1=0.
	r2=Step
	Integral=0.
	while r2 <= r:
		Value=air((r2+r1)/2., Wavelength, Mirror_diameter)
		Integral=Integral+Value*pi*(r2*r2-r1*r1)
		r1 = r2
		r2 = r2 + Step
	return (Integral)

def zyl(Step, Limit):
	r1=0.
	r2=Step
	Integral=0.
	while r2 <= Limit:
		#Value=pow((r2+r1)/2.,2.)
		Value=1
		#Integral=Integral+Value*Parameter*Arcsec_multiplier*Parameter*Arcsec_multiplier*pi*(r2*r2-r1*r1)
		#print Parameter*Arcsec_multiplier*Parameter*Arcsec_multiplier*pi*(r2*r2-r1*r1)
		Integral=Integral+Value*pi*(r2*r2-r1*r1)
		r1 = r2
		r2 = r2 + Step
	return (Integral)
	

def moffat(r, Beta, Seeing):
	Alpha=Seeing/(2.*sqrt(pow(2,(1/Beta))-1))
	moffat=(Beta-1)/(pi*Alpha*Alpha)*pow((1+(r*r)/(Alpha*Alpha)), (-Beta))
	return(moffat)

def moffat_int(r, Beta, Seeing):
	parameter_b=Beta
	parameter_a=Seeing/(2.*sqrt(pow(2,(1/parameter_b))-1))
	parameter_c=r
	integral=(1.-1./(pow(1.+(parameter_c*parameter_c)/(parameter_a*parameter_a), parameter_b-1.)))
	return(integral)

def moffat_func(x, y, Beta, Seeing):
	Alpha=Seeing/(2.*sqrt(pow(2.,(1./Beta))-1.))
	moffat=(Beta-1.)/(pi*Alpha*Alpha)*pow((1.+(x*x+y*y)/(Alpha*Alpha)), (-Beta))
	return(moffat)

def moffat_func_int(Beta, Limit_x1, Limit_x2, Limit_y1, Limit_y2, Seeing):
	Step=abs(Limit_x1-Limit_x2)/50.
	x=Limit_x1
	y=Limit_y1
	Integral=0.
	while x<Limit_x2:
		while y<Limit_y2:
			Value=moffat_func(x, y, Beta, Seeing)
			Integral=Integral+Value*Step*Step
			#Temp="%s\t%s\t%s\n" % (x, y, Funktion)
			#File.write(Temp)
			y=y+Step
		y=Limit_y1
		x=x+Step
	return(Integral)
	
def gauss_func(x,y,Sigma):
	Temp = 1./(2.0*math.pi*Sigma**2)*exp(-0.5*(x**2+y**2)/Sigma**2)
	return(Temp)

def gauss_func_int(Limit_x1, Limit_x2, Limit_y1, Limit_y2, Seeing):
	Step=abs(Limit_x1-Limit_x2)/50.
	x=Limit_x1
	y=Limit_y1
	Integral=0.
	while x<Limit_x2:
		while y<Limit_y2:
			Sigma = Seeing / (2*sqrt(2*log(2)))
			Value=gauss_func(x, y, Sigma)
			Integral=Integral+Value*Step*Step
			y=y+Step
		y=Limit_y1
		x=x+Step
	return(Integral)

def gauss_func_int_lowres(Limit_x1, Limit_x2, Limit_y1, Limit_y2, Seeing):
	Step=abs(Limit_x1-Limit_x2)/10.
	x=Limit_x1
	y=Limit_y1
	Integral=0.
	while x<Limit_x2:
		while y<Limit_y2:
			Sigma = Seeing / (2.*sqrt(2.*log(2.)))
			Value=gauss_func(x, y, Sigma)
			Integral=Integral+Value*Step*Step
			y=y+Step
		y=Limit_y1
		x=x+Step
	return(Integral)

def gauss_int(Radius,Seeing):
	Temp = 0.0
	Temp1 = 0.0
	Step = Radius/100.
	Sigma = Seeing / (2*sqrt(2*log(2)))
	x = -Radius
	y = -Radius
	while x <= Radius:
		#print x
		while y <= Radius:
			if math.sqrt(x**2 + y**2) < Radius:
				Temp = gauss_func(x,y,Sigma)*Step*Step
				Temp1 = Temp1 + Temp
				Temp = Temp1
			y = y + Step
		y = -Radius
		x = x + Step
	return(Temp)


def psf_func_int(Beta, Limit_x1, Limit_x2, Limit_y1, Limit_y2, Seeing, Wavelength, Mirror_diameter, Strehl):
	x=Limit_x1
	y=Limit_y1
	Integral=0.
	Norm_moffat = moffat_func(0.0,0.0, Beta, Seeing)
	Norm_air = air_func(0.0,0.0,Wavelength, Mirror_diameter)
	F0 =  (Norm_moffat-Strehl*Norm_air) / (Norm_moffat - Norm_air)
	Integral = (1.-F0)*moffat_func_int(Beta, Limit_x1, Limit_x2, Limit_y1, Limit_y2, Seeing)
	Integral = Integral + F0*air_func_int(Limit_x1, Limit_x2, Limit_y1, Limit_y2, Wavelength, Mirror_diameter)
	#print "xy",F0
	#print "xy-moffat", (1.-F0)*moffat_func_int(Beta, Limit_x1, Limit_x2, Limit_y1, Limit_y2, Seeing)
	#print "xy-air", F0*air_func_int(Limit_x1, Limit_x2, Limit_y1, Limit_y2, Wavelength, Mirror_diameter)
	#print F0*air_func_int(Limit_x1, Limit_x2, Limit_y1, Limit_y2, Wavelength, Mirror_diameter)
	#File=open('test.csv','w')
	#while x<Limit_x2:
		#while y<Limit_y2:
			#Part_1=moffat_func(x, y, Beta, Seeing)
			#Integral=Integral+Part_1*Step*Step
			#Part_2=air_func(x, y, Wavelength, Mirror_diameter)
			#Integral=Integral+((Beta-1.)/(pi*Alpha*Alpha*Strehl))*Part_2*Step*Step
			##Temp="%.9e\t%.9e\t%.9e\t%.9e\n" % (x, y, Part_1, Part_2)
			##File.write(Temp)
			##print Integral
			##Temp="%.9e\t%.9e\t%.9e\n" % (x, y, Teil_2)
			##Temp="%.9e\t%.9e\t%.9e\n" % (x*360.*3600./(2.*pi), y*360.*3600./(2.*pi), Teil_2)
			#y=y+Step
		#y=Limit_y1
		#x=x+Step
	#File.close()
	return(Integral)

def psf_int(r, Beta, Seeing, Wavelength, Mirror_diameter, Strehl):
	Norm_moffat = moffat(0.0, Beta, Seeing)
	Norm_air = air(0.0,Wavelength, Mirror_diameter)
	F0 =  (Norm_moffat-Strehl*Norm_air) / (Norm_moffat - Norm_air)
	Integral=(1.-F0)*moffat_int(r, Beta, Seeing) + F0*air_int(r, Wavelength, Mirror_diameter)
	#print "r",F0
	#print "r-moffat",(1.-F0)*moffat_int(r, Beta, Seeing)
	#print "r-air",F0*air_int(r, Wavelength, Mirror_diameter)
	return(Integral)

#def slit_int(Wavelength, Slit_width, limit1, limit2):
#	slit_function=pow(sin(pi/Wavelength*Slit_width*sin(x))/(Wavelength*Slit_width*sin(x)),2)
#	return(slit_width)

def linear_interpolation(Filename, Datacontainer_def, Datacontainer, Wavelength_step_new, Redshift):
#	Filename_w='test.csv'
	Table_1={}
	Table_1_def={}
	String_1=""
	String_2=""
	Fileobject_read = open(Filename, 'r')
#	Fileobject_write = open(Filename_w, 'w')
	Data = Fileobject_read.readline()
	Row = string.split(Data)
	Number_of_datalines = int(Row[0])
	Wavelength_start = (float(Row[1]))*Redshift + (float(Row[1]))
#	Wavelength_start = (float(Row[1]))
	Wavelength_end = (float(Row[2]))*Redshift + (float(Row[2]))
#	Wavelength_end = (float(Row[2]))
	Wavelength_step = (1+Redshift)*(float(Row[3]))
#	Wavelength_step = (float(Row[3]))
	if Wavelength_step_new==0:
		Wavelength_step_new=Wavelength_step
	Table_1[Filename]=[Filename, Number_of_datalines, Wavelength_start, Wavelength_end, Wavelength_step]
	Wavelength_start_new = int((Wavelength_start-.2)/Wavelength_step_new)*Wavelength_step_new+0.2
	if Wavelength_start_new<Wavelength_start:
		Wavelength_start_new=Wavelength_start_new+Wavelength_step_new
	Wavelength_end_new = int((Wavelength_end-.2)/Wavelength_step_new)*Wavelength_step_new+0.2+Wavelength_step_new
	if Wavelength_end_new>Wavelength_end:
		Wavelength_end_new=Wavelength_end_new-Wavelength_step_new
	String = "%.10f" % ((Wavelength_end_new-Wavelength_start_new)/Wavelength_step_new)
	Number_of_datalines_new = int(float(String))+1
	Datacontainer[0]=Filename
	Datacontainer[1]=Number_of_datalines_new
	Datacontainer[2]=Wavelength_start_new
	Datacontainer[3]=Wavelength_end_new
	Datacontainer[4]=Wavelength_step_new
	String ="%d\t%.5f\t%.5f\t%.5f\n" % (Number_of_datalines_new, Wavelength_start_new, Wavelength_end_new, Wavelength_step_new)
#	Fileobject_write.write(String)	
	for i in range(Number_of_datalines):
		Data = Fileobject_read.readline()
		Row = string.split(Data)
		String="%.9f" % (float(Row[0])*Redshift + float(Row[0]))
		Table_1_def[String] = float(Row[1])
	for i in range(Number_of_datalines_new):
		String="%.9f" % (i*Wavelength_step_new+Wavelength_start_new)
		Datacontainer_def[String]=0.0
	Temp=0.0
	k=0
	i=0
	if Wavelength_step_new > Wavelength_step:
		for j in range(1,Number_of_datalines_new):
			String_j="%.9f" % (j*Wavelength_step_new+Wavelength_start_new)
			Wavelength_j = j*Wavelength_step_new+Wavelength_start_new
			Wavelength_i = Wavelength_start+i*Wavelength_step
			String_i = "%.9f" % (Wavelength_i)
			while Wavelength_i < Wavelength_j-Wavelength_step_new/2.:
				i=i+1
				Wavelength_i=Wavelength_i+Wavelength_step
			String_i = "%.9f" % (Wavelength_i)
			while Wavelength_i >= Wavelength_j-Wavelength_step_new/2. and Wavelength_i <= Wavelength_j+Wavelength_step_new/2. and Wavelength_i <= Wavelength_end:
				Temp=Temp+Table_1_def[String_i]
				i=i+1
				k=k+1
				Wavelength_i=Wavelength_i+Wavelength_step
				String_i = "%.9f" % (Wavelength_i)
			if k>0:
				Datacontainer_def[String_j]=Temp/k
			i=i-1
			Temp=0.0
			k=0
	elif Wavelength_step_new <= Wavelength_step:
		j=0
		for i in range(1,Number_of_datalines):
			Wavelength_i = (i*Wavelength_step+Wavelength_start)
			Wavelength_i_minus_1 = ((i-1)*Wavelength_step+Wavelength_start)
			String_i = "%.9f" % (i*Wavelength_step+Wavelength_start)
			String_i_minus_1 = "%.9f" % Wavelength_i_minus_1
			Wavelength_j = j*Wavelength_step_new+Wavelength_start_new
			String_j="%.9f" % Wavelength_j
			while j <= Number_of_datalines_new and i < Number_of_datalines and Wavelength_j <= Wavelength_i:
				M=(Table_1_def[String_i]-Table_1_def[String_i_minus_1])/(Wavelength_i-Wavelength_i_minus_1)
				B=(Table_1_def[String_i]-M*Wavelength_i)
				Datacontainer_def[String_j]=M*Wavelength_j+B
				j=j+1
				Wavelength_i = (i*Wavelength_step+Wavelength_start)
				Wavelength_i_minus_1 = ((i-1)*Wavelength_step+Wavelength_start)
				String_i = "%.9f" % Wavelength_i
				String_i_minus_1 = "%.9f" % Wavelength_i_minus_1
				Wavelength_j = j*Wavelength_step_new+Wavelength_start_new
				String_j="%.9f" % Wavelength_j
	Fileobject_read.close()
#	for i in range(Number_of_datalines_new):
#		String_j="%.9f" % (i*Wavelength_step_new+Wavelength_start_new)
#		String="%.9f\t%.9f\n" % (float(String_j), Datacontainer_def[String_j])
#		Fileobject_write.write(String)
#	Fileobject_write.close()
	return(Datacontainer_def, Datacontainer)

def linear_interpolation_output(Filename, Datacontainer_def, Datacontainer, Wavelength_step_new, Redshift):
	######################
	Filename_w='test_alt.csv'
	######################
	Table_1={}
	Table_1_def={}
	String_1=""
	String_2=""
	Fileobject_read = open(Filename, 'r')
	######################
	Fileobject_write = open(Filename_w, 'w')
	######################
	Data = Fileobject_read.readline()
	Row = string.split(Data)
	Number_of_datalines = int(Row[0])
	Wavelength_start = (float(Row[1]))*Redshift + (float(Row[1]))
#	Wavelength_start = (float(Row[1]))
	Wavelength_end = (float(Row[2]))*Redshift + (float(Row[2]))
#	Wavelength_end = (float(Row[2]))
	Wavelength_step = (1+Redshift)*(float(Row[3]))
#	Wavelength_step = (float(Row[3]))
	if Wavelength_step_new==0:
		Wavelength_step_new=Wavelength_step
	Table_1[Filename]=[Filename, Number_of_datalines, Wavelength_start, Wavelength_end, Wavelength_step]
	Wavelength_start_new = int((Wavelength_start-.2)/Wavelength_step_new)*Wavelength_step_new+0.2
	if Wavelength_start_new<Wavelength_start:
		Wavelength_start_new=Wavelength_start_new+Wavelength_step_new
	Wavelength_end_new = int((Wavelength_end-.2)/Wavelength_step_new)*Wavelength_step_new+0.2+Wavelength_step_new
	if Wavelength_end_new>Wavelength_end:
		Wavelength_end_new=Wavelength_end_new-Wavelength_step_new
	String = "%.10f" % ((Wavelength_end_new-Wavelength_start_new)/Wavelength_step_new)
	Number_of_datalines_new = int(float(String))+1
	Datacontainer[0]=Filename
	Datacontainer[1]=Number_of_datalines_new
	Datacontainer[2]=Wavelength_start_new
	Datacontainer[3]=Wavelength_end_new
	Datacontainer[4]=Wavelength_step_new
	String ="%d\t%.10f\t%.10f\t%.10f\n" % (Number_of_datalines_new, Wavelength_start_new, Wavelength_end_new, Wavelength_step_new)
	######################
	Fileobject_write.write(String)	
	######################
	for i in range(Number_of_datalines):
		Data = Fileobject_read.readline()
		Row = string.split(Data)
		String="%.10f" % (float(Row[0])*Redshift + float(Row[0]))
		Table_1_def[String] = float(Row[1])
	for i in range(Number_of_datalines_new):
		String="%.10f" % (i*Wavelength_step_new+Wavelength_start_new)
		Datacontainer_def[String]=0.0
	Temp=0.0
	k=0
	i=0
	if Wavelength_step_new > Wavelength_step:
		for j in range(1,Number_of_datalines_new):
			String_j="%.10f" % (j*Wavelength_step_new+Wavelength_start_new)
			Wavelength_j = j*Wavelength_step_new+Wavelength_start_new
			Wavelength_i = Wavelength_start+i*Wavelength_step
			String_i = "%.10f" % (Wavelength_i)
			while Wavelength_i < Wavelength_j-Wavelength_step_new/2.:
				i=i+1
				Wavelength_i=Wavelength_i+Wavelength_step
			String_i = "%.10f" % (Wavelength_i)
			while Wavelength_i >= Wavelength_j-Wavelength_step_new/2. and Wavelength_i <= Wavelength_j+Wavelength_step_new/2. and Wavelength_i <= Wavelength_end:
				Temp=Temp+Table_1_def[String_i]
				i=i+1
				k=k+1
				Wavelength_i=Wavelength_i+Wavelength_step
				String_i = "%.10f" % (Wavelength_i)
			if k>0:
				Datacontainer_def[String_j]=Temp/k
			i=i-1
			Temp=0.0
			k=0
	elif Wavelength_step_new <= Wavelength_step:
		j=0
		for i in range(1,Number_of_datalines):
			Wavelength_i = (i*Wavelength_step+Wavelength_start)
			Wavelength_i_minus_1 = ((i-1)*Wavelength_step+Wavelength_start)
			String_i = "%.10f" % (i*Wavelength_step+Wavelength_start)
			String_i_minus_1 = "%.10f" % Wavelength_i_minus_1
			Wavelength_j = j*Wavelength_step_new+Wavelength_start_new
			String_j="%.10f" % Wavelength_j
			while j <= Number_of_datalines_new and i < Number_of_datalines and Wavelength_j <= Wavelength_i:
				M=(Table_1_def[String_i]-Table_1_def[String_i_minus_1])/(Wavelength_i-Wavelength_i_minus_1)
				B=(Table_1_def[String_i]-M*Wavelength_i)
				Datacontainer_def[String_j]=M*Wavelength_j+B
				j=j+1
				Wavelength_i = (i*Wavelength_step+Wavelength_start)
				Wavelength_i_minus_1 = ((i-1)*Wavelength_step+Wavelength_start)
				String_i = "%.10f" % Wavelength_i
				String_i_minus_1 = "%.10f" % Wavelength_i_minus_1
				Wavelength_j = j*Wavelength_step_new+Wavelength_start_new
				String_j="%.10f" % Wavelength_j
	Fileobject_read.close()
	######################
	for i in range(Number_of_datalines_new):
		String_j="%.10f" % (i*Wavelength_step_new+Wavelength_start_new)
		String="%.10f\t%.10f\n" % (float(String_j), Datacontainer_def[String_j])
		Fileobject_write.write(String)
	Fileobject_write.close()
	######################
	return(Datacontainer_def, Datacontainer)

def background_theoretical(Atmosphere_def, Airmass, Filename, Datacontainer_def, Wavelength_start, Wavelength_end, Wavelength_step):
	#Wavelength_step=0.001
	#print Airmass
	Temp=[]
	List=[]
	String=[]
	Limit=[]
	Factor=[]
	Sum=[]
	OH=0.0
	#Atmosphere_def={}
	#Atmosphere=['air_Atmosphere/16_150.csv', 0, 0.0, 0.0, 0.0]
	#List.sort()
	#print List[0]
	Atmosphere_temperature=250.
	Zodiacal_temperature=5800.
	#Wavelength_start=0.801
	#Wavelength_start = int((0.8-.2)/Wavelength_step)*Wavelength_step+0.2
	#if Wavelength_start<0.8:
	#	Wavelength_start=Wavelength_start+Wavelength_step
	#Wavelength_end = int((2.5-.2)/Wavelength_step)*Wavelength_step+0.2+Wavelength_step
	#if Wavelength_end>2.5:
	#	Wavelength_end=Wavelength_end-Wavelength_step
	h=6.63e-34
	c=299830000.
	k=1.3807e-23
	Intensity_tmp=0.0
	Wavelength=0.0
	Wavelength=Wavelength_start
	Number_of_datalines = int((Wavelength_end-Wavelength_start)/Wavelength_step+1)
	ini_read = open('background.ini', 'r')
	Temp = ini_read.readlines()
	#print Temp
	for i in range(1, len(Temp)):
		Raw = string.split(Temp[i])
		#print Raw
		Limit.append(float(Raw[1]))
		Factor.append(float(Raw[2]))
		Sum.append(float(Raw[3]))
#	Fileobject_write = open('background.csv', 'w')
	Database = open(Filename, 'r')
	Data=Database.readline()
	Number_of_lines = len(Database.readlines())
	Database.close()
	Database = open(Filename, 'r')
	#linear_interpolation('air_transmission/16_150.csv', Atmosphere_def, Atmosphere, Wavelength_step)
	#Liste=Atmosphere_def.keys()
	#Liste.sort()
#	for i in Liste:
#		print i
	Data=Database.readline()
	Data=Database.readline()
	Row=string.split(Data)
	OH=float(Row[0])
#	Temp="%d\t%.9f\t%.9f\t%.9f\n" % (Number_of_datalines, Wavelength_start, Wavelength_end, Wavelength_step)
	Datacontainer=[Filename, Number_of_datalines, Wavelength_start, Wavelength_end, Wavelength_step]
#	Fileobject_write.write(Temp)
	i=1
	j=0
	while i:
		if (Wavelength-Wavelength_end) > 1e-9:
			i=0
			break
		while (OH < Wavelength+Wavelength_step) and (OH >= Wavelength):
			if Wavelength < Limit[j]:
				#print Wavelength, Limit[j], Factor[j], Sum[j]
				Intensity_tmp=Intensity_tmp+float(Row[1])*Factor[j]+Sum[j]
				#Intensity_tmp=Intensity_tmp+float(Row[1])
				Data = Database.readline()
				if not Data :
					break
				Row = string.split(Data)
				OH=float(Row[0])
			else:
				j=j+1
		if (OH < Wavelength) and Data:
			Data = Database.readline()
			if Data :
				Row = string.split(Data)
				OH=float(Row[0])
		else:	
			String="%.9f" % (Wavelength)
			Blackbody_atmosphere = 4.*h*c*c/pow((Wavelength/1000000.),5)/(exp(h*c/((Wavelength/1000000.)*k*Atmosphere_temperature))-1.)
			Blackbody_Zodiacal = 120.*2.*h*c*c/pow((Wavelength/1000000.),5)/(exp(h*c/((Wavelength/1000000.)*k*Zodiacal_temperature))-1.)/5e14
			Intensity_tmp = Intensity_tmp*Atmosphere_def[String]+(1.-Atmosphere_def[String])*Blackbody_atmosphere+Blackbody_Zodiacal
			Intensity_tmp = Intensity_tmp*(-0.000278719*Airmass*Airmass*Airmass-0.0653841*Airmass*Airmass+1.11979*Airmass-0.0552132)
			#Temp="%.10f\t%e\t%.10f\n" % (Wavelength, Intensity_tmp, Atmosphere_def[String])
#			Temp="%.9f\t%.9f\t%.4e\n" % (Wavelength, Intensity_tmp, Atmosphere_def[String])
#			Fileobject_write.write(Temp)
			Temp="%.9f" % (Wavelength)
			Datacontainer_def[Temp]=Intensity_tmp
			Intensity_tmp=0.0
			Wavelength=Wavelength+Wavelength_step
	Database.close()
#	Fileobject_write.close()
	return (Datacontainer_def, Datacontainer)

def background_lines(Filename, Datacontainer_def, Wavelength_start, Wavelength_end, Wavelength_step):
	Temp=[]
	List=[]
	String=[]
	Limit=[]
	Factor=[]
	Sum=[]
	OH=0.0
	Intensity_tmp=0.0
	Wavelength=0.0
	Wavelength=Wavelength_start
	Number_of_datalines = int((Wavelength_end-Wavelength_start)/Wavelength_step+1)
	ini_read = open('background.ini', 'r')
	Temp = ini_read.readlines()
	for i in range(1, len(Temp)):
		Raw = string.split(Temp[i])
		Limit.append(float(Raw[1]))
		Factor.append(float(Raw[2]))
		Sum.append(float(Raw[3]))
	Database = open(Filename, 'r')
	Data=Database.readline()
	Number_of_lines = len(Database.readlines())
	Database.close()
	Database = open(Filename, 'r')
	Data=Database.readline()
	Data=Database.readline()
	Row=string.split(Data)
	OH=float(Row[0])
	Datacontainer=[Filename, Number_of_datalines, Wavelength_start, Wavelength_end, Wavelength_step]
	i=1
	j=0
	while i:
		if (Wavelength-Wavelength_end) > 1e-9:
			i=0
			break
		while (OH < Wavelength+Wavelength_step) and (OH >= Wavelength):
			if Wavelength < Limit[j]:
				Intensity_tmp=Intensity_tmp+float(Row[1])*Factor[j]+Sum[j]
				Data = Database.readline()
				if not Data :
					break
				Row = string.split(Data)
				OH=float(Row[0])
			else:
				j=j+1
		if (OH < Wavelength) and Data:
			Data = Database.readline()
			if Data :
				Row = string.split(Data)
				OH=float(Row[0])
		else:	
			String="%.9f" % (Wavelength)
			Datacontainer_def[String]=Intensity_tmp
			Intensity_tmp=0.0
			Wavelength=Wavelength+Wavelength_step
	Database.close()
	return (Datacontainer_def, Datacontainer)


def background_theoretical_rev(Atmosphere_def, Airmass, Datacontainer_def, Wavelength_start, Wavelength_end, Wavelength_step):
	Temp=[]
	List=[]
	String=[]
	Limit=[]
	Factor=[]
	Sum=[]
	OH=0.0
	Atmosphere_temperature=250.
	Zodiacal_temperature=5800.
	h=6.63e-34
	c=299830000.
	k=1.3807e-23
	Intensity_tmp=0.0
	Wavelength=Wavelength_start
	Number_of_datalines = int((Wavelength_end-Wavelength_start)/Wavelength_step+1)
	for i in range(0,len(Datacontainer_def)):
		String="%.9f" % (Wavelength)
		Blackbody_atmosphere = 4.*h*c*c/pow((Wavelength/1000000.),5)/(exp(h*c/((Wavelength/1000000.)*k*Atmosphere_temperature))-1.)
		Blackbody_Zodiacal = 120.*2.*h*c*c/pow((Wavelength/1000000.),5)/(exp(h*c/((Wavelength/1000000.)*k*Zodiacal_temperature))-1.)/5e14
		Datacontainer_def[String] = Datacontainer_def[String]*Atmosphere_def[String]+(1.-Atmosphere_def[String])*Blackbody_atmosphere+Blackbody_Zodiacal
		Datacontainer_def[String] = Datacontainer_def[String]*(-0.000278719*Airmass*Airmass*Airmass-0.0653841*Airmass*Airmass+1.11979*Airmass-0.0552132)
		Temp="%.9f" % (Wavelength)
		Wavelength=Wavelength+Wavelength_step
	return (Datacontainer_def)