procedure dbs (r_images,b_images,comb_list,r_bias,b_bias,r_dflat,b_dflat,
	over_fn,t_order,trim,trim_sec,bias_sec,wid_no,
	b_fun,b_order,b_sampl,tr_fn,tr_order,arcs,flux_arcs,
        flux_std,f_comb_list,smooth,sm_arcs,sm_f_arcs,
	sm_std,sm_comb_list,sm_f_std,sm_f_comb_list,
	cal_dir,red_data,red_comb)

string r_images {prompt="File of red images to be reduced"}
string b_images {prompt="File of blue images to be reduced"}
string comb_list {prompt="Combination file for images"}
string r_bias	{prompt="File of red bias frames to be combined"}
string b_bias   {prompt="File of blue bias frames to be combined"}
string r_dflat	{prompt="File of red dome flat frames to be combined"}
string b_dflat	{prompt="File of blue dome flat frames to be combined"}
string over_fn {"legendre",prompt="Function to fit trace of overscan region"}
int    t_order {2,prompt="Order of overscan fitting function"}
string bias_sec {"[1793:1838,1:200]",prompt="Overscan region to fit"}
bool trim    {yes,prompt="Trim image?"}
string trim_sec {"[30:1720,1:200]",prompt="Trim section"}
int wid_no {1,prompt="Number of different widths to extract"}
string b_fun {"legendre",prompt="Background function"}
int b_order {1,prompt="Background function order"}
string b_sampl {"-80:-20,20:80",prompt="Background sample regions"}
string tr_fn {"spline3",prompt="Trace fitting function"}
int tr_order {3,prompt="Trace fitting function order"}
string arcs {prompt="File of object arcs"}
string flux_arcs {prompt="File of object flux standard arcs"}
string flux_std {prompt="File of flux standards"}
string f_comb_list {prompt="Combination file for flux standards"}
bool smooth {yes,prompt="Sky absorption correct?"}
string sm_arcs {prompt="File of object smooth spectrum standard arcs"}
string sm_f_arcs {prompt="File of flux smooth spectrum standard arcs"}
string sm_std {prompt="File of object smooth spectrum standards"}
string sm_comb_list {prompt="Combination file for object smooth spectrum standards"}
string sm_f_std {prompt="File of flux smooth spectrum standards"}
string sm_f_comb_list {prompt="Combination file for flux smooth spectrum standards"}
string cal_dir {"/priv/madras4/lkewley/cal/",prompt="directory containing flux calibration data"}
string red_data {prompt="has data already been bias subtracted and flatfielded?"}
string red_comb {prompt="has data already been combined?"}

struct *ld_r_images
struct *ld_b_images
real *wid_ld_1
int *ld_comb_list
int *ld_f_comb_list
int *ld_sm_comb_list
string *arcs_ld
string *flux_arcs_ld
string *sm_arcs_ld
string *flux_std_ld
string  *sm_std_ld
string *rimage_c_ld
string *bimage_c_ld
string *image_c_ld
string *r_flux_ld
string *b_flux_ld
string *stand_ld
string *r_sm_ld
string *b_sm_ld
string *r_text_ld
string *b_text_ld
string *sm_c_ld
string *flux_c_ld
string *r_sm_temp_ld
string *b_sm_temp_ld
string *r_sm_list_ld
string *b_sm_list_ld	
string *r_coefs_ld
string *b_coefs_ld
string *images_c_ld
string *cal_name_ld
string *sm_f_std_ld
string *sm_f_arcs_ld

begin
	string r_images_l,rimage,bimage,rc_name,bc_name,rflux,bflux
	string b_images_l,bias_sec_l,red,aper_l,blue,r_flux,cname
	string r_bias_l,comb_lis_l,b_bias_l,r_arc_out,b_arc_out
	string r_dflat_l,b_dflat_l,arc_list_l,r_ref,b_ref,bc_name_o
	string over_fn_l,no2test,r_arc_name,b_arc_name,rc_name_o
	string trim_sec_l,b_fun_l,b_sampl_l,tr_fn_l,blue_comb
	string r_sm,b_sm,b_flux,rarc,barc,rsm,bsm,imageclist,rcomb
	string redfluxlist,bluefluxlist,frc_name,fbc_name,frc_name_o
	string fbc_name_o,rsmlist,bsmlist,sm_rc_name,sm_bc_name,bcomb
	string sm_rc_name_o,sm_bc_name_o,fluxclist,smclistfile,rsmtemplist
	string bsmtemplist,r_text,b_text,r_coefs,b_coefs,cal_name
	string rclist,bclist
	int t_order_l,tr_order_l,aper_no_l,b_order_l,clist,a,b,c,obj_no 
	int fclist,smclist,flux_no,sm_no,flux_flag,sm_flag,comb_sm_flag
	int comb_flux_flag,total_no,flag,abs_flag,flag2,sens_flag,sm_f_no
	int airmass_flag
	bool trim_l
	real a1,a2,a3,half_a1,half_a2,half_a3,a1_min,r_airmass[10]
	real a2_min,a3_min,b_airmass[10],r_ave_airmass,b_ave_airmass
	real r_airmass_o,b_airmass_o,r_coef1,r_coef2,b_coef1,b_coef2
	real r_airmass_no,b_airmass_no,rlambda,blambda
	real r_exposed,b_exposed,r_obj_exposed,b_obj_exposed
	file obj_list_blue,obj_list_red,rflux.list,bflux.list,rarc.list
	file barc.list,red_o1.list,blue_o1.list,red_o2.list,blue_o2.list
	file red_o3.list,blue_o3.list,rimage_c.list,bimage_c.list
	file image_c.list,rflux_c.list,bflux_c.list,flux_c.list,rsm_c.list
	file bsm_c.list,sm_c.list,stand.list,r_flux.list,b_flux.list 
	file r_sm.list,b_sm.list,b_sm_temp.list,r_sm_temp.list
	file rsm_t.list,bsm_t.list,sm_c_t.list,airmass.txt
			
	# define list-directed parameters
	ld_r_images = r_images
	ld_b_images = b_images
	ld_comb_list=comb_list
	arcs_ld=arcs
	flux_arcs_ld=flux_arcs
	sm_arcs_ld=sm_arcs
	flux_std_ld=flux_std
	sm_std_ld=sm_std
	ld_f_comb_list=f_comb_list
	ld_sm_comb_list=sm_comb_list
	sm_f_std_ld=sm_f_std
	sm_f_arcs_ld=sm_f_arcs

# check that the necessary packages are loaded
	if (! defpac("noao") || ! defpac ("onedspec") || ! defpac ("twodspec") || ! defpac("apextract") || ! defpac("imred") || ! defpac ("ccdred"))
	bye
	
# check existance of input files
	if (!access(r_images)){
		print((r_images),' not found, aborting')
		beep ; beep ; goto cleanup}
	if (!access(b_images)){
		print((b_images),' not found, aborting')
		beep ; beep ; goto cleanup}
	if ((red_comb)=='no'){
	if (!access(comb_list)){
		print((comb_list),' not found, aborting')
		beep ; beep ; goto cleanup}
	}
	if ((red_data)=='no'){
	if (!access(r_bias)){
		print((r_bias),' not found, aborting')
		beep ; beep ; goto cleanup}
	if (!access(b_bias)){
		print((b_bias),' not found, aborting')
		beep ; beep ; goto cleanup}
	if (!access(r_dflat)){
		print((r_dflat),' not found, aborting')
		beep ; beep ; goto cleanup}
	if (!access(b_dflat)){
		print((b_dflat),' not found, aborting')
		beep ; beep ; goto cleanup}
	}
	if (!access(arcs)){
		print((arcs),' not found, aborting')
		beep ; beep ; goto cleanup}
	if (!access(sm_arcs)){
		print((sm_arcs),' not found, aborting')
		beep ; beep ; goto cleanup}
	if (!access(flux_arcs)){
		print((flux_arcs),' not found, aborting')
		beep ; beep ; goto cleanup}
	if (!access(sm_std)){
		print((sm_std),' not found, aborting')
		beep ; beep ; goto cleanup}
	if ((red_comb)=='no'){
	if (!access(sm_comb_list)){
		print((sm_comb_list),' not found, aborting')
		beep ; beep ; goto cleanup}
	}	
	if (!access(flux_std)){
		print((flux_std),' not found, aborting')
		beep ; beep ; goto cleanup}
	if ((red_comb)=='no'){
	if (!access(f_comb_list)){
		print((f_comb_list),' not found, aborting')
		beep ; beep ; goto cleanup}
	}

# noao.twodspec.apextract.dispaxis=1	

# get the query parameters
	r_images_l=r_images
	b_images_l=b_images
	r_bias_l = r_bias
	b_bias_l= b_bias
	r_dflat_l=r_dflat
	b_dflat_l=b_dflat
	over_fn_l=over_fn
	t_order_l=t_order
	bias_sec_l=bias_sec
	trim_l=trim
	trim_sec_l=trim_sec	
	b_fun_l=b_fun
	b_order_l=b_order
	b_sampl_l=b_sampl
	tr_fn_l = tr_fn
	tr_order_l = tr_order

if ((red_data)=='yes'){
	goto red_data_stage1
}

#overscan correct bias frames
print('overscan correcting bias frames')
ccdproc('@'//r_bias_l,fixp-,overs+,trim-,zeroc-,darkc-,flatc-,
	illumc-,fringec-,readc-,scanc-,readax="line",biass=(bias_sec_l),
	intera+,functio=(over_fn_l),order=(t_order_l),ccdtype="",
        )
ccdproc('@'//b_bias_l,fixp-,overs+,trim-,zeroc-,darkc-,flatc-,
	illumc-,fringec-,readc-,scanc-,readax="line",biass=(bias_sec_l),
	intera+,functio=(over_fn_l),order=(t_order_l),ccdtype="",
        )


# combine bias frames
 print('combining bias frames')
if (access('R_bias_c.fits')){
	print('red bias frames have already been combined')
	goto step2}
 imcombine ('@'//r_bias_l, outp="R_bias_c",combine="median", logf="STDOUT", 
	rej="avsigclip", proj-,outt="real",offs="none",maskty="none",
	scale="none",zero="none",weight="none",stat="[100:1700,10:190]",
	mclip+)
step2:
if (access('B_bias_c.fits')){
	print('blue bias frames have already been combined')
	goto step3}
 imcombine ('@'//b_bias_l, outp="B_bias_c", combine="median", logf="STDOUT",
	rej="avsigclip", proj-,outt="real",offs="none",maskty="none",
	scale="none",zero="none",weight="none",stat="[100:1700,10:190]",
	mclip+)
step3:

# overscan correcting and bias subtraction of dome flats
print('overscan correcting and bias subtraction of dome flats')
ccdproc('@'//r_dflat_l,fixp-,overs+,trim=(trim_l),zeroc+,darkc-,flatc-,
	illumc-,fringec-,readc-,scanc-,readax="line",biass=(bias_sec_l),
	trims=(trim_sec_l),zero="R_bias_c",intera+,
	functio=(over_fn_l),order=(t_order_l),ccdtype="")
	
ccdproc('@'//b_dflat_l,fixp-,overs+,trim=(trim_l),zeroc+,darkc-,flatc-,
	illumc-,fringec-,readc-,scanc-,readax="line",biass=(bias_sec_l),
	trims=(trim_sec_l),zero="B_bias_c",intera+,
	functio=(over_fn_l),order=(t_order_l),ccdtype="")


# combine dome flats
print('combining dome flats')
if (access('R_dflat_c.fits')){
	print('red dome flats have already been combined')
	goto step4}
 imcombine ('@'//r_dflat_l,outp="R_dflat_c",combine="median",logf="STDOUT",
	rej="avsigclip",proj-,outt="real",offs="none",maskty="none",
	scale="mode",zero="none",weight="none",stat="[100:1700,10:190]",
	mclip+)
step4:	
if (access('B_dflat_c.fits')){
	print('blue dome flats have already been combined')
	goto step5}
imcombine ('@'//b_dflat_l, outp="B_dflat_c", combine="median", logf="STDOUT",	rej="avsigclip", proj-,outt="real",offs="none",maskty="none",
	scale="mode",zero="none",weight="none",stat="[100:1700,10:190]",
	mclip+)

# apply bias correction, flatfield correction, overscan strip 
# correction, and trim data.

step5:
print('Bias, flatfield, and overscan correction for red images')

 ccdproc('@'//r_images_l,fixp-,overs+,trim=(trim_l),zeroc+,darkc-,flatc+,
	illumc-,fringec-,readc-,scanc-,readax="line",biass=(bias_sec_l),
	trims=(trim_sec_l),zero="R_bias_c",flat="R_dflat_c",intera+,
	functio=(over_fn_l),order=(t_order_l),ccdtype="")

print('Bias, flatfield, and overscan correction for blue images')

 ccdproc('@'//b_images_l,fixp-,overs+,trim=(trim_l),zeroc+,darkc-,flatc+,
	illumc-,fringec-,readc-,scanc-,readax="line",biass=(bias_sec_l),
	trims=(trim_sec_l),zero="B_bias_c",flat="B_dflat_c",intera+,
	functio=(over_fn_l),order=(t_order_l),ccdtype="")

# apply bias correction, flatfield correction, overscan strip correction
# and trim standards

red_data_stage1:

if (access('rflux.list')){
		delete rflux.list
}
if (access('bflux.list')){
		delete bflux.list
}

flux_std_ld=(flux_std)

while (fscan (flux_std_ld,rflux,bflux,cal_name) !=EOF){
	print ((rflux), >> 'rflux.list')
	print ((bflux), >> 'bflux.list')
}

if (access('rsm.list')){
		delete rsm.list
}
if (access('bsm.list')){
		delete bsm.list
}
if (access('rsm_t.list')){
		delete rsm_t.list
}
if (access('bsm_t.list')){
		delete bsm_t.list
}
if (access('rsm_f.list')){
		delete rsm_f.list
}
if (access('bsm_f.list')){
		delete bsm_f.list
}
while (fscan (sm_std_ld,rsm,bsm) !=EOF){
	print ((rsm), >> 'rsm.list')
	print ((bsm), >> 'bsm.list')
	print ((rsm), >> 'rsm_t.list')
	print ((bsm), >> 'bsm_t.list')
}


sm_f_std_ld=(sm_f_std)
while (fscan (sm_f_std_ld,rsm,bsm) !=EOF){
	print ((rsm), >> 'rsm_t.list')
	print ((bsm), >> 'bsm_t.list')
	print ((rsm), >> 'rsm_f.list')
	print ((bsm), >> 'bsm_f.list')
}

if ((red_data)=='yes'){
	goto red_data_stage2
}

print('Bias, flatfield, and overscan correction for red smooth spectrum standards')

ccdproc('@rsm_t.list',fixp-,overs+,trim=(trim_l),zeroc+,darkc-,flatc+,
	illumc-,fringec-,readc-,scanc-,readax="line",biass=(bias_sec_l),
	trims=(trim_sec_l),zero="R_bias_c",flat="R_dflat_c",intera+,
	functio=(over_fn_l),order=(t_order_l),ccdtype="")

print('Bias, flatfield, and overscan correction for blue smooth spectrum standards')

ccdproc('@bsm_t.list',fixp-,overs+,trim=(trim_l),zeroc+,darkc-,flatc+,
	illumc-,fringec-,readc-,scanc-,readax="line",biass=(bias_sec_l),
	trims=(trim_sec_l),zero="B_bias_c",flat="B_dflat_c",intera+,
	functio=(over_fn_l),order=(t_order_l),ccdtype="")


print('Bias, flatfield, and overscan correction for red flux standards')

ccdproc('@rflux.list',fixp-,overs+,trim=(trim_l),zeroc+,darkc-,flatc+,
	illumc-,fringec-,readc-,scanc-,readax="line",biass=(bias_sec_l),
	trims=(trim_sec_l),zero="R_bias_c",flat="R_dflat_c",intera+,
	functio=(over_fn_l),order=(t_order_l),ccdtype="")

print('Bias, flatfield, and overscan correction for blue flux standards')

ccdproc('@bflux.list',fixp-,overs+,trim=(trim_l),zeroc+,darkc-,flatc+,
	illumc-,fringec-,readc-,scanc-,readax="line",biass=(bias_sec_l),
	trims=(trim_sec_l),zero="B_bias_c",flat="B_dflat_c",intera+,
	functio=(over_fn_l),order=(t_order_l),ccdtype="")

red_data_stage2:

if (access('rarc.list')){
		delete rarc.list
}
if (access('barc.list')){
		delete barc.list
}
if (access('airmass.txt')){
		delete 'airmass.txt'
}
while (fscan (arcs_ld,rarc,barc) !=EOF){
	print ((rarc), >> 'rarc.list')
	print ((barc), >> 'barc.list')
}

if ((red_data)=='yes'){
	goto red_data_stage3
}

print('Bias, flatfield, and overscan correction for red arcs')

ccdproc('@rarc.list',fixp-,overs+,trim=(trim_l),zeroc+,darkc-,flatc+,
	illumc-,fringec-,readc-,scanc-,readax="line",biass=(bias_sec_l),
	trims=(trim_sec_l),zero="R_bias_c",flat="R_dflat_c",intera+,
	functio=(over_fn_l),order=(t_order_l),ccdtype="")

print('Bias, flatfield, and overscan correction for blue arcs')

ccdproc('@barc.list',fixp-,overs+,trim=(trim_l),zeroc+,darkc-,flatc+,
	illumc-,fringec-,readc-,scanc-,readax="line",biass=(bias_sec_l),
	trims=(trim_sec_l),zero="B_bias_c",flat="B_dflat_c",intera+,
	functio=(over_fn_l),order=(t_order_l),ccdtype="")

red_data_stage3:

if ((red_comb)=='yes'){
	goto red_data_stage4
}

a=0

if (access('image_c.list')){
		delete 'image_c.list'
}
if (access('rimage_c.list')){
		delete 'rimage_c.list'
}
if (access('bimage_c.list')){
		delete 'bimage_c.list'
}

obj_no=0
ld_r_images=r_images
ld_b_images=b_images
airmass_flag=0

while (fscan(ld_comb_list,clist)!=EOF){
	
	a=a+1
	
	b=0
	c=0
	if (access ('red_ctemp.list')){
		delete 'red_ctemp.list'
	}
	if (access ('blue_ctemp.list')){
		delete 'blue_ctemp.list'
	}

	for (i=1; i <= 10; i+=1){
		r_airmass[i]=0.
		b_airmass[i]=0.
	}

	for (i=1; i <= (clist); i+=1){
		rimage=ld_r_images
		bimage=ld_b_images 


		print ((rimage), >> 'red_ctemp.list')
		print ((bimage), >> 'blue_ctemp.list')	

		if (i == 1) {
			obj_no=obj_no+1
			rc_name=(rimage//'_C')
			bc_name=(bimage//'_C')
			rc_name_o=(rimage)
			bc_name_o=(bimage)
			print((rimage),' ',(bimage), >> 'image_c.list')
			print((rimage), >> 'rimage_c.list')
			print((bimage), >> 'bimage_c.list')	
			
		}		

		if (access(rc_name//'.fits')){
			goto step6}

#		get airmass and calculate average for each red object

		imgets(rimage,param='airmass')
		r_airmass[i]=real(imgets.value)
		imgets(bimage,param='airmass')
		b_airmass[i]=real(imgets.value)
		
		print((rimage),' ',(r_airmass[i]),' ',(bimage),' ',(b_airmass[i]),>>'airmass.txt')	
	
step6:
	
	}
	if (access(rc_name//'.fits')){
		goto step7}

	r_airmass_o=0.
	r_ave_airmass=0.
	b_airmass_o=0.
	b_ave_airmass=0.
	r_airmass_no=0	
	b_airmass_no=0


	for (i=1; i<=10; i+=1){
		if (r_airmass[i] != 0.){
		if (r_airmass[i] <3.){
			r_airmass_no=r_airmass_no+1
			r_ave_airmass=r_airmass[i]+r_airmass_o
 			r_airmass_o=r_ave_airmass		
		}
		}
	}
	if ((r_airmass_no)!=0.){
	r_ave_airmass=r_airmass_o/r_airmass_no
	}
	if ((r_airmass_no)==0.){
	r_ave_airmass=0.
	}

	for (i=1; i<=10; i+=1){
		if (b_airmass[i] != 0.){
		if (b_airmass[i] <3.){
			b_airmass_no=b_airmass_no+1
			b_ave_airmass=b_airmass[i]+b_airmass_o
 			b_airmass_o=b_ave_airmass
		}	
		}
	}
	if ((b_airmass_no)!=0.){
	b_ave_airmass=b_airmass_o/b_airmass_no
	}
	if ((b_airmass_no)==0.) {
	b_ave_airmass=0.
	}
	
	print('average red airmass= ',(r_ave_airmass),'average blue airmass= ',(b_ave_airmass),>>'airmass.txt')	
	
	if (r_ave_airmass==0.){
		r_ave_airmass=b_ave_airmass
	}
	if (b_ave_airmass==0.){
		b_ave_airmass=r_ave_airmass
	}
	if (r_ave_airmass==0.){
		airmass_flag=1
	}
	if (b_ave_airmass==0.){
		airmass_flag=1
	}
	if (r_ave_airmass > 3.0){
		airmass_flag=1
	}
	if (b_ave_airmass > 3.0) {
		airmass_flag=1
	}
	
	
step7:
	if (clist != 1) {	

	if (access(rc_name//'.fits')){
	print('red images for object ',(a),' have already been combined')
	goto step8}
	
	print("Combining images with CR rejection for Object ",(a))
	imcombine (inp='@red_ctemp.list',output=rc_name,
	logf="STDOUT",combine="median",proje-,outty="real",
	reject="avsigclip",scale="mode",mclip+,stat="[100:700,10:190]")	
	hedit(images=rc_name,fields='airmass',value=r_ave_airmass,
	add+,del-,verify-,show-,update+)
step8:
	if (access(bc_name//'.fits')){
	print('blue images for object ',(a),' have already been combined')
	goto step9}
	imcombine (inp='@blue_ctemp.list',output=bc_name,
	logf="STDOUT",combine="median",proje-,outty="real",
	reject="avsigclip",scale="mode",mclip+,stat="[100:700,10:190]")	
	hedit(images=bc_name,fields='airmass',value=b_ave_airmass,
	add+,del-,verify-,show-,update+)	
step9:
	}
	else{
	if (access(rc_name//'.fits')){
	print('red images for object ',(a),' have already been combined')
	goto imcomb_step10}
	imcopy(input=rc_name_o,output=rc_name,verbose-)
imcomb_step10:
	if (access(bc_name//'.fits')){
	print('blue images for object ',(a),' have already been combined')
	goto imcomb_step11}
	imcopy(input=bc_name_o,output=bc_name,verbose-)
imcomb_step11:
	}
}




# Combining flux standards

if (access('r_flux.list')){
		delete r_flux.list
}
if (access('b_flux.list')){
		delete b_flux.list
}
if (access('cal_name.list')){
		delete cal_name.list
}

flux_std_ld=flux_std
while (fscan(flux_std_ld,rflux,bflux,cal_name)!=EOF){
	print((rflux), >> 'r_flux.list')
	print((bflux), >> 'b_flux.list')
	print((cal_name), >> 'cal_name.list')
}
	
a=0

redfluxlist='r_flux.list'
bluefluxlist='b_flux.list'
r_flux_ld=redfluxlist
b_flux_ld=bluefluxlist
cal_name_ld='cal_name.list'
flux_no=0	

if (access('flux_c.list')){
		delete 'flux_c.list'
}
if (access('rflux_c.list')){
		delete 'rflux_c.list'
}
if (access('bflux_c.list')){
		delete 'bflux_c.list'
}

rflux=''
bflux=''
ld_f_comb_list=f_comb_list

while (fscan(ld_f_comb_list,fclist)!=EOF){
	a=a+1
	b=0
	c=0
	
	if(access('frcomb.list')){
		delete 'frcomb.list'
	}	
	if(access('fbcomb.list')){
		delete 'fbcomb.list'
	}	
	for (i=1; i<=10; i+=1){
		r_airmass[i]=0.
		b_airmass[i]=0.
	}
	for (i=1; i<=fclist; i+=1){
		rflux=r_flux_ld
		bflux=b_flux_ld
		cname=cal_name_ld

		print((rflux), >> 'frcomb.list')
		print((bflux), >> 'fbcomb.list')
		

		if (i == 1) {
			flux_no=flux_no+1
			frc_name=(rflux//'_C')
			fbc_name=(bflux//'_C')
			frc_name_o=(rflux)
			fbc_name_o=(bflux)
			print((rflux),' ',(bflux),' ',(cname),>> 'flux_c.list')
			print((rflux), >> 'rflux_c.list')
			print((bflux), >> 'bflux_c.list')	

		}
		if (access((frc_name)//'.fits')){
		print((rflux),' has already been combined, skipping to next standard')
				goto fcomb_step1
			}
#		get airmass and calculate average for each object

		imgets(rflux,param='airmass')
		r_airmass[i]=real(imgets.value)
		imgets(bflux,param='airmass')
		b_airmass[i]=real(imgets.value)
		
		print((rflux),' ',(r_airmass[i]),' ',(bflux),' ',(b_airmass[i]),>>'airmass.txt')	

				
fcomb_step1:	
	}
	if (access((frc_name)//'.fits')){
				goto fcomb_step2}

	r_airmass_o=0.
	r_ave_airmass=0.
	b_airmass_o=0.
	b_ave_airmass=0.
	r_airmass_no=0	
	b_airmass_no=0

	for (i=1; i<=10; i+=1){
		if (r_airmass[i] != 0.){
		if (r_airmass[i] <3. ){
			r_airmass_no=r_airmass_no+1
			r_ave_airmass=r_airmass[i]+r_airmass_o
 			r_airmass_o=r_ave_airmass		
		}
		}
	}
	if ((r_airmass_no)!=0.){
	r_ave_airmass=r_airmass_o/r_airmass_no
	}
	if ((r_airmass_no)==0.){
	r_ave_airmass=0.
	}

	for (i=1; i<=10; i+=1){
		if (b_airmass[i] != 0.){
		if (b_airmass[i] <3.){
			b_airmass_no=b_airmass_no+1
			b_ave_airmass=b_airmass[i]+b_airmass_o
 			b_airmass_o=b_ave_airmass
		}	
		}
	}
	
	if ((b_airmass_no)!=0){
	b_ave_airmass=b_airmass_o/b_airmass_no
	}
	if (b_airmass_no==0){
	b_ave_airmass=0.
	}
	
	
	print((rimage),' ',(r_ave_airmass),' ',(bimage),' ',(b_ave_airmass),>>'airmass.txt')	
	
	if (r_ave_airmass==0.){
		r_ave_airmass=b_ave_airmass
	}
	if (b_ave_airmass==0.){
		b_ave_airmass=r_ave_airmass
	}
	if (r_ave_airmass==0.){
		airmass_flag=1
	}
	if (b_ave_airmass==0.){
		airmass_flag=1
	}
	if (r_ave_airmass > 3.0){
		airmass_flag=1
	}
	if (b_ave_airmass > 3.0) {
		airmass_flag=1
	}
	

fcomb_step2:	

	if (fclist != 1) {	

	if (access((frc_name)//'.fits')){
				goto fcomb_step3}
	print("Combining images with CR rejection for Flux Standard ",(a))
		
			imcombine (inp='@frcomb.list',output=frc_name,
			logf="STDOUT",combine="median",proje-,outty="real",
			reject="avsigclip",scale="mode",mclip+,
			stat="[100:700,10:190]")	
			hedit(images=frc_name,fields='airmass',
			value=r_ave_airmass,add+,del-,verify-,show-,update+)
fcomb_step3:
			if (access((fbc_name)//'.fits')){
		print((bflux),' has already been combined, skipping to next standard')
				goto fcomb_step4
			}
			imcombine (inp='@fbcomb.list',output=fbc_name,
			logf="STDOUT",combine="median",proje-,outty="real",
			reject="avsigclip",scale="mode",mclip+,
			stat="[100:700,10:190]")	
			hedit(images=fbc_name,fields='airmass',
			value=b_ave_airmass,add+,del-,verify-,show-,update+)
fcomb_step4:		
		}
	if (fclist == 1){
		if (access((frc_name)//'.fits')){
		print((rflux),' has already been combined, skipping to next standard')
				goto fcomb_step5
			}
		imcopy(input=frc_name_o,output=frc_name,verbose-)
fcomb_step5:
		if (access((fbc_name)//'.fits')){
		print((bflux),' has already been combined, skipping to next standard')
				goto fcomb_step6
			}
		imcopy(input=fbc_name_o,output=fbc_name,verbose-)
fcomb_step6:
	}

flux_comb_step1:

}

if (airmass_flag==1){
print('odd airmass, see airmass.txt for information, exiting')
goto cleanup
}


#  Combining smooth spectrum standards


rsmlist='rsm.list'
bsmlist='bsm.list'
r_sm_ld=rsmlist
b_sm_ld=bsmlist


a=0
sm_no=0	
if (access('sm_c.list')){
		delete 'sm_c.list'
}
if (access('sm_c_t.list')){
		delete 'sm_c_t.list'
}
if (access('rsm_c.list')){
		delete 'rsm_c.list'
}
if (access('bsm_c.list')){		
	delete 'bsm_c.list'
}

rsm=''
bsm=''


while (fscan(ld_sm_comb_list,smclist)!=EOF){
	a=a+1
	b=0
	c=0
	if(access('sm_rcomb.list')){
		delete 'sm_rcomb.list'
	}	
	if(access('sm_bcomb.list')){
		delete 'sm_bcomb.list'
	}	

	for (i=1; i<=smclist; i+=1){
		
		rsm=r_sm_ld
		bsm=b_sm_ld
		print((rsm), >> 'sm_rcomb.list')
		print((bsm), >> 'sm_bcomb.list')

		if (i == 1) {
			sm_no=sm_no+1
			sm_rc_name=(rsm//'_C')
			sm_bc_name=(bsm//'_C')
			sm_rc_name_o=(rsm)
			sm_bc_name_o=(bsm)
			print((rsm),' ',(bsm), >> 'sm_c.list')
			print((rsm), >> 'rsm_c.list')
			print((bsm), >> 'bsm_c.list')
			print((rsm),' ',(bsm), >> 'sm_c_t.list')	
						
		}		
	}
	
	print("Combining images with CR rejection for Smooth Spectrum Standard ",(a))
		
		if (access((sm_rc_name)//'.fits')){
	print('Smooth spectrum standard ',(a),' has already been combined, skipping to next standard')
			goto sm_comb_step1
 		}

	if (smclist != 1) {			
			imcombine (inp='@sm_rcomb.list',output=sm_rc_name,
			logf="STDOUT",combine="median",proje-,outty="real",
			reject="avsigclip",scale="mode",mclip+,
			stat="[100:700,10:190]")	

			imcombine (inp='@sm_bcomb.list',output=sm_bc_name,
			logf="STDOUT",combine="median",proje-,outty="real",
			reject="avsigclip",scale="mode",mclip+,
			stat="[100:700,10:190]")	
		
	}
	if (smclist == 1){
		imcopy(input=sm_rc_name_o,output=sm_rc_name,verbose-)
		imcopy(input=sm_bc_name_o,output=sm_bc_name,verbose-)
	}

sm_comb_step1:
}

rsmlist='rsm_f.list'
bsmlist='bsm_f.list'
r_sm_ld=rsmlist
b_sm_ld=bsmlist


a=0
sm_f_no=0	
if (access('sm_c_f.list')){
		delete 'sm_c_f.list'
}
if (access('rsm_c_f.list')){
		delete 'rsm_c_f.list'
}
if (access('bsm_c_f.list')){		
	delete 'bsm_c_f.list'
}

rsm=''
bsm=''

ld_sm_comb_list=(sm_f_comb_list)

while (fscan(ld_sm_comb_list,smclist)!=EOF){
	a=a+1
	b=0
	c=0
	if(access('sm_rcomb.list')){
		delete 'sm_rcomb.list'
	}	
	if(access('sm_bcomb.list')){
		delete 'sm_bcomb.list'
	}	

	for (i=1; i<=smclist; i+=1){
		
		rsm=r_sm_ld
		bsm=b_sm_ld
		print((rsm), >> 'sm_rcomb.list')
		print((bsm), >> 'sm_bcomb.list')

		if (i == 1) {
			sm_f_no=sm_f_no+1
			sm_rc_name=(rsm//'_C')
			sm_bc_name=(bsm//'_C')
			sm_rc_name_o=(rsm)
			sm_bc_name_o=(bsm)
			print((rsm),' ',(bsm), >> 'sm_c_f.list')
			print((rsm), >> 'rsm_c_f.list')
			print((bsm), >> 'bsm_c_f.list')	
			print((rsm),' ',(bsm), >> 'sm_c_t.list')			
		}		
	}

	print("Combining images with CR rejection for Smooth Spectrum Standard ",(a))
		
		if (access((sm_rc_name)//'.fits')){
	print('Smooth spectrum standard ',(a),' has already been combined, skipping to next standard')
			goto sm_comb_step2
 		}
	if (smclist != 1) {	
			imcombine (inp='@sm_rcomb.list',output=sm_rc_name,
			logf="STDOUT",combine="median",proje-,outty="real",
			reject="avsigclip",scale="mode",mclip+,
			stat="[100:700,10:190]")	

			imcombine (inp='@sm_bcomb.list',output=sm_bc_name,
			logf="STDOUT",combine="median",proje-,outty="real",
			reject="avsigclip",scale="mode",mclip+,
			stat="[100:700,10:190]")	
		
	}
	if (smclist ==1){
	imcopy(input=sm_rc_name_o,output=sm_rc_name,verbose-)
	imcopy(input=sm_bc_name_o,output=sm_bc_name,verbose-)
	}
sm_comb_step2:
}

red_data_stage4:

if ((red_comb)=='yes'){

	ld_r_images=(r_images)
	ld_b_images=(b_images)
	sm_std_ld=(sm_std)
	sm_f_std_ld=(sm_f_std)
	flux_std_ld=(flux_std)

	if(access('image_c.list')){
		delete 'image_c.list'
	}
	if(access('sm_c_f.list')){
		delete 'sm_c_f.list'
	}
	if(access('rsm_c_f.list')){
		delete 'rsm_c_f.list'
	}
	if(access('bsm_c_f.list')){
		delete 'bsm_c_f.list'
	}	
	if(access('sm_c_t.list')){
		delete 'sm_c_t.list'
	}
	if(access('sm_c.list')){
		delete 'sm_c.list'
	}
	if(access('rsm_c.list')){
		delete 'rsm_c.list'
	}
	if(access('bsm_c.list')){
		delete 'bsm_c.list'
	}
	if(access('sm_c_t.list')){
		delete 'sm_c_t.list'
	}
	if(access('flux_c.list')){
		delete 'flux_c.list'
	}
	if(access('rflux_c.list')){
		delete 'rflux_c.list'
	}
	if(access('bflux_c.list')){
		delete 'bflux_c.list'
	}
	obj_no=0
	sm_no=0
	sm_f_no=0
	flux_no=0
	while(fscan (ld_r_images, red)!=EOF && fscan(ld_b_images, blue)!=EOF){
		print((red),' ',(blue), >> 'image_c.list')
		obj_no=obj_no+1
	}
	while(fscan (sm_f_std_ld, rsm, bsm)!=EOF){
			sm_f_no=sm_f_no+1
			print((rsm),' ',(bsm), >> 'sm_c_f.list')
			print((rsm), >> 'rsm_c_f.list')
			print((bsm), >> 'bsm_c_f.list')	
			print((rsm),' ',(bsm), >> 'sm_c_t.list')	
	}

	while(fscan (sm_std_ld, rsm, bsm)!=EOF){
			sm_no=sm_no+1
			print((rsm),' ',(bsm), >> 'sm_c.list')
			print((rsm), >> 'rsm_c.list')
			print((bsm), >> 'bsm_c.list')
			print((rsm),' ',(bsm), >> 'sm_c_t.list')	
	}			
	while( fscan(flux_std_ld, rflux,bflux,cname)!=EOF){ 
	print((rflux),' ',(bflux),' ',(cname),>> 'flux_c.list')
			flux_no=flux_no+1
			print((rflux), >> 'rflux_c.list')
			print((bflux), >> 'bflux_c.list')	
	}
}


imageclist='image_c.list'
image_c_ld=imageclist
ld_r_images=(r_images)
ld_b_images=(b_images)
sm_std_ld=(sm_std)
sm_f_std_ld=(sm_f_std)
flux_std_ld=(flux_std)

# image extraction 

	a=0
 while (fscan (image_c_ld, red, blue) !=EOF){

	a=a+1
	
	for(j=1; j<=wid_no; j+=1){

	print('Extracting Red images for object ',(a),' aperture ',(j))

	if (access(red//'_o'//j//'.fits')){
	print('Object ',(red),' aperture ',(j),' has already been extracted')
		goto step10
	}
	apall(red//'_C',outp=red//'_o'//j,inter+,find+,recen+,resiz-,edit+,trace+,
		fittra+,extrac+,extras+,review+,b_func=(b_fun_l),
		b_order=(b_order_l),b_sampl=(b_sampl_l),nsum=50,nfind=1,
		t_func=(tr_fn_l),t_order=(tr_order_l),back="median",
		skyb=1,weight="none",pfit="fit1d",readno=0,
		gain=1,lsigma=4,usigma=4,nsubaps=1)
step10:

	print('Extracting Blue images for object ',(a),' aperture ',(j))

	if (access(blue//'_o'//j//'.fits')){
	print ('Object ',(blue),' aperture ',(j),' has already been extracted')
	goto step11
	}

	apall(blue//'_C',outp=blue//"_o"//j,inter+,find+,recen+,resiz-,edit+,trace+,
		fittra+,extrac+,extras+,review+,b_func=(b_fun_l),
		b_order=(b_order_l),b_sampl=(b_sampl_l),nsum=50,nfind=1,
		t_func=(tr_fn_l),t_order=(tr_order_l),back="median",
		skyb=1,weight="none",pfit="fit1d",readno=0,
		gain=1,lsigma=4,usigma=4,nsubaps=1)
step11:
	}
}

# Extraction of flux and smooth spectrum standards

print('Extracting standard stars')
	
if (access('std.list')){
	delete 'std.list'
}

sm_c_ld='sm_c_t.list'
fluxclist='flux_c.list'
flux_c_ld=fluxclist
r_sm=''
b_sm=''


while (fscan (sm_c_ld,r_sm,b_sm)!=EOF && fscan(flux_c_ld,r_flux,b_flux,cal_name)!=EOF){
	print((r_flux),' ',(b_flux),' ',(r_sm),' ',(b_sm), >> 'std.list')
}
	
smclistfile='sm_c_t.list'
sm_c_ld=smclistfile
fluxclist='flux_c.list'
flux_c_ld=fluxclist
r_sm=''
b_sm=''

while (fscan (flux_c_ld,r_flux,b_flux,cal_name)!=EOF){
	if (access(r_flux//"_o.fits")){
		print('already extracted ',r_flux//"_o",' skipping to next standard')
		goto sm_extraction
	}
		print('extracting ',(r_flux))
		apall(r_flux//'_C',outp=r_flux//"_o",inter+,find+,recen+,resiz+,edit+,trace+,
		fittra+,extrac+,extras+,review+,b_func=(b_fun_l),
		b_order=(b_order_l),b_sampl=(b_sampl_l),nfind=1,
		t_func=(tr_fn_l),t_order=(tr_order_l),back="median",
		skyb=1,weight="none",pfit="fit1d",readno=0,nfind=1,
		gain=1,lsigma=4,usigma=4,nsubaps=1,nsum=10)

	print('extracting ',(b_flux))
	apall(b_flux//'_C',outp=b_flux//"_o",inter+,find+,recen+,resiz+,edit+,trace+,
		fittra+,extrac+,extras+,review+,b_func=(b_fun_l),
		b_order=(b_order_l),b_sampl=(b_sampl_l),nfind=1,
		t_func=(tr_fn_l),t_order=(tr_order_l),back="median",
		skyb=1,weight="none",pfit="fit1d",readno=0,nfind=1,
		gain=1,lsigma=4,usigma=4,nsubaps=1,nsum=10)	

sm_extraction:	
}

while (fscan (sm_c_ld,r_sm,b_sm)!=EOF){
	
	if (access(r_sm//"_o.fits")){
		print('already extracted ',r_sm//"_o",' skipping to next standard')
		goto sm_step1
	}
	

	print('extracting ',(r_sm))
	apall(r_sm//'_C',outp=r_sm//"_o",inter+,find+,recen+,resiz+,edit+,trace+,
		fittra+,extrac+,extras+,review+,b_func=(b_fun_l),
		b_order=(b_order_l),b_sampl=(b_sampl_l),nfind=1,
		t_func=(tr_fn_l),t_order=(tr_order_l),back="median",
		skyb=1,weight="none",pfit="fit1d",readno=0,nfind=1,
		gain=1,lsigma=4,usigma=4,nsubaps=1,nsum=10)
sm_step1:
if (access(b_sm//"_o.fits")){
		print('already extracted ',b_sm//"_o",' skipping to next standard')
		goto sm_step2
	}
	print('extracting ',(b_sm))
	apall(b_sm//'_C',outp=b_sm//"_o",inter+,find+,recen+,resiz+,edit+,trace+,
		fittra+,extrac+,extras+,review+,b_func=(b_fun_l),
		b_order=(b_order_l),b_sampl=(b_sampl_l),nfind=1,
		t_func=(tr_fn_l),t_order=(tr_order_l),back="median",
		skyb=1,weight="none",pfit="fit1d",readno=0,nfind=1,
		gain=1,lsigma=4,usigma=4,nsubaps=1,nsum=10)
sm_step2:	
}	


	red=''
	blue=''

#  Wavelength Calibration	

	arcs_ld=arcs
	imageclist='image_c.list'
	image_c_ld=imageclist

while (fscan(image_c_ld, red, blue)!=EOF && fscan(arcs_ld,r_arc_name,b_arc_name,r_ref,b_ref)!=EOF){	

	if (access(red//'_w'//1//'.fits')){
	print('Object ',(red),' has already been wavelength calibrated')
	goto step12}
	b=0
	for (i=1 ; i<=obj_no ; i+=1){
		if (access(r_arc_name//"_t"//b//'.fits')){
			b=b+1
		}
	}
	r_arc_out=r_arc_name//"_t"//b
	b_arc_out=b_arc_name//"_t"//b
	rc_name=red//'_o'
	bc_name=blue//'_o'

	
	print("Wavelength Calibrating Object ",(red)," and ",(blue))


	apall(r_arc_name,out=(r_arc_out)//'.fits',ref=(red)//'_C',recen-,trace-,
	back-,intera-,t_fun='spline3')
	
	apall(b_arc_name,out=(b_arc_out)//'.fits',ref=(blue)//'_C',recen-,trace-,
	back-,intera-,t_fun='spline3')
		

	if ((r_arc_out) == (r_ref//"_t0")){
		print ('identifying ',(r_arc_name)) 


		identify(images=r_arc_out,secti='middle line',datab='database',
		coordli='linelists$idhenear.dat',nsum=10,match=10.,
		maxfeat=75,zwidth=100.,ftype='emission',fwid=2.,cradius=4.,
		thres=20.,mins=2.,functi='legendre',order=2,niter=0,low_rej=3.,
		high_rej=3.,grow=0.,autow-,graph='stdgraph',autowri+)
		print ('identifying ',(b_arc_name)) 

		identify(images=b_arc_out,secti='middle line',datab='database',
		coordli='linelists$idhenear.dat',nsum=10,match=10.,
		maxfeat=75,zwidth=100.,ftype='emission',fwid=2.,cradius=4.,
		thres=20.,mins=2.,functi='legendre',order=2,niter=0,low_rej=3.,
		high_rej=3.,grow=0.,autow-,graph='stdgraph',autowri+)
		

	}

	if ((r_arc_out) != (r_ref//'_t0')){
	print ('reidentifying ',(r_arc_name), ' using reference ',(r_ref))
 
		reidentify(images=r_arc_out,refer=r_ref//'_t0',interac+,
		secti='middle line',newaps+,overr+,refit+,trace-,step=10,
		nsum=10,shift=0.,nlost=0,crad=5.,thresh=10.,add-,
		coord='linelists$idhenear.dat',match=10.,maxf=50,mins=2.,
		datab='database',logf='logfile',verb-,graph='stdgraph',ans+)

	print ('reidentifying ',(b_arc_name), ' using reference ',(b_ref))
 
		reidentify(images=b_arc_out,refer=b_ref//'_t0',interac+,
		secti='middle line',newaps+,overr+,refit+,trace-,step=10,
		nsum=10,shift=0.,nlost=0,crad=5.,thresh=10.,add-,
		coord='linelists$idhenear.dat',match=10.,maxf=50,mins=2.,
		datab='database',logf='logfile',verb-,graph='stdgraph',ans+)
	
	}

	for (i=1; i<=wid_no;i+=1){

		hedit(images=rc_name//i,fields='REFSPEC1',
		value=r_arc_out,add+,del-,ver-,show-,update+)
		hedit(images=bc_name//i,fields='REFSPEC1',
		value=b_arc_out,add+,del-,ver-,show-,update+)
		dispcor(inp=rc_name//i,outp=red//'_w'//i,lin+,log-,flux+,
		same-,glob-,ignor-,confirm+,listo-,verb+,datab='database')
		dispcor(inp=bc_name//i,outp=blue//'_w'//i,lin+,log-,flux+,
		same-,glob-,ignor-,confirm+,listo-,verb+,datab='database')

	}	
step12:
}


#  Wavelength Calibration of Smooth Spectrum Standards	

sm_c_ld='sm_c.list'
sm_arcs_ld=sm_arcs

while (fscan(sm_c_ld, red, blue)!=EOF &&
fscan(sm_arcs_ld,r_arc_name,b_arc_name,r_ref,b_ref)!=EOF){	
	
	
	total_no=obj_no+sm_no
	b=0
	for (i=1 ; i<=total_no ; i+=1){
		if (access(r_arc_name//"_t"//b//'.fits')){
			b=b+1
		}
	}
	r_arc_out=r_arc_name//"_t"//b
	b_arc_out=b_arc_name//"_t"//b
	rc_name=red//'_o'
	bc_name=blue//'_o'

	
	print("Wavelength Calibrating Smooth Spectrum Standard ",(red)," and ",(blue))

	if (access(red//'_w.fits')){
		print((red),' has already been wavelength calibrated, skipping to next standard')
		goto sm_w_step1
		
	 }
		apall(r_arc_name,out=(r_arc_out)//'.fits',ref=(red)//'_C',
		recen-,trace-,back-,intera-,t_fun='spline3')
		if ((r_arc_out) == (r_ref//'_t0')){
		print ('identifying ',(r_arc_name)) 
		identify(images=r_arc_out,secti='middle line',datab='database',
		coordli='linelists$idhenear.dat',nsum=10,match=10.,
		maxfeat=75,zwidth=100.,ftype='emission',fwid=2.,cradius=4.,
		thres=20.,mins=2.,functi='legendre',order=2,niter=0,low_rej=3.,
		high_rej=3.,grow=0.,autow-,graph='stdgraph',autowri+)
		}
	if ((r_arc_out) != (r_ref//'_t0')){
	print ('reidentifying ',(r_arc_name), ' using reference ',(r_ref))
 
		reidentify(images=r_arc_out,refer=r_ref//'_t0',interac+,
		secti='middle line',newaps+,overr+,refit+,trace-,step=10,
		nsum=10,shift=0.,nlost=0,crad=5.,thresh=10.,add-,
		coord='linelists$idhenear.dat',match=10.,maxf=50,mins=2.,
		datab='database',logf='logfile',verb-,graph='stdgraph',ans+)
	}
	hedit(images=rc_name,fields='REFSPEC1',
	value=r_arc_out,add+,del-,ver-,show-,update+)
	dispcor(inp=rc_name,outp=red//'_w',lin+,log-,flux+,
	same-,glob-,ignor-,confirm+,listo-,verb+)
		
	
sm_w_step1:
		if (access(blue//'_w.fits')){
		print((blue),' has already been wavelength calibrated, skipping to next standard')
		goto sm_w_step2
		
	 }
		apall(b_arc_name,out=(b_arc_out)//'.fits',ref=(blue)//'_C',
		recen-,trace-,back-,intera-,t_fun='spline3')
		

	if ((b_arc_out) == (b_ref//'_t0')){
		
		print ('identifying ',(b_arc_name)) 
		identify(images=b_arc_out,secti='middle line',datab='database',
		coordli='linelists$idhenear.dat',nsum=10,match=10.,
		maxfeat=75,zwidth=100.,ftype='emission',fwid=2.,cradius=4.,
		thres=20.,mins=2.,functi='legendre',order=2,niter=0,low_rej=3.,
		high_rej=3.,grow=0.,autow-,graph='stdgraph',autowri+)
		

	}

	

	if ((b_arc_out) != (b_ref//'_t0')){
	
	print ('reidentifying ',(b_arc_name), ' using reference ',(b_ref))
 
		reidentify(images=b_arc_out,refer=b_ref//'_t0',interac+,
		secti='middle line',newaps+,overr+,refit+,trace-,step=10,
		nsum=10,shift=0.,nlost=0,crad=5.,thresh=10.,add-,
		coord='linelists$idhenear.dat',match=10.,maxf=50,mins=2.,
		datab='database',logf='logfile',verb-,graph='stdgraph',ans+)
	
	}

		
		hedit(images=bc_name,fields='REFSPEC1',
		value=b_arc_out,add+,del-,ver-,show-,update+)
		dispcor(inp=bc_name,outp=blue//'_w',lin+,log-,flux+,
		same-,glob-,ignor-,confirm+,listo-,verb+)

sm_w_step2:	
}

sm_c_ld='sm_c_f.list'
sm_arcs_ld=(sm_f_arcs)

while (fscan(sm_c_ld, red, blue)!=EOF && fscan(sm_arcs_ld,r_arc_name,b_arc_name,r_ref,b_ref)!=EOF){	
	sm_flag=0
	if (access(red//'_w.fits')){
		print((red),' has already been wavelength calibrated, skipping to next standard')
		sm_flag=1
		
	 }
	
	if (sm_flag==0){
		total_no=obj_no+sm_no+sm_f_no
		b=0
		for (i=1 ; i<=total_no ; i+=1){
			if (access(r_arc_name//"_t"//b//'.fits')){
				b=b+1
			}
		}
		r_arc_out=r_arc_name//"_t"//b
		b_arc_out=b_arc_name//"_t"//b
		rc_name=red//'_o'
		bc_name=blue//'_o'

	
	print("Wavelength Calibrating Smooth Spectrum Standard ",(red)," and ",(blue))

		apall(r_arc_name,out=(r_arc_out)//'.fits',ref=(red)//'_C',
		recen-,trace-,back-,intera-,t_fun='spline3')
	
		apall(b_arc_name,out=(b_arc_out)//'.fits',ref=(blue)//'_C',
		recen-,trace-,back-,intera-,t_fun='spline3')
		

	if ((r_arc_out) == (r_ref//'_t0')){
		print ('identifying ',(r_arc_name)) 
		identify(images=r_arc_out,secti='middle line',datab='database',
		coordli='linelists$idhenear.dat',nsum=10,match=10.,
		maxfeat=75,zwidth=100.,ftype='emission',fwid=2.,cradius=4.,
		thres=20.,mins=2.,functi='legendre',order=2,niter=0,low_rej=3.,
		high_rej=3.,grow=0.,autow-,graph='stdgraph',autowri+)

		print ('identifying ',(b_arc_name)) 
		identify(images=b_arc_out,secti='middle line',datab='database',
		coordli='linelists$idhenear.dat',nsum=10,match=10.,
		maxfeat=75,zwidth=100.,ftype='emission',fwid=2.,cradius=4.,
		thres=20.,mins=2.,functi='legendre',order=2,niter=0,low_rej=3.,
		high_rej=3.,grow=0.,autow-,graph='stdgraph',autowri+)
		

	}

	if ((r_arc_out) != (r_ref//'_t0')){
	print ('reidentifying ',(r_arc_name), ' using reference ',(r_ref))
 
		reidentify(images=r_arc_out,refer=r_ref//'_t0',interac+,
		secti='middle line',newaps+,overr+,refit+,trace-,step=10,
		nsum=10,shift=0.,nlost=0,crad=5.,thresh=10.,add-,
		coord='linelists$idhenear.dat',match=10.,maxf=50,mins=2.,
		datab='database',logf='logfile',verb-,graph='stdgraph',ans+)

	print ('reidentifying ',(b_arc_name), ' using reference ',(b_ref))
 
		reidentify(images=b_arc_out,refer=b_ref//'_t0',interac+,
		secti='middle line',newaps+,overr+,refit+,trace-,step=10,
		nsum=10,shift=0.,nlost=0,crad=5.,thresh=10.,add-,
		coord='linelists$idhenear.dat',match=10.,maxf=50,mins=2.,
		datab='database',logf='logfile',verb-,graph='stdgraph',ans+)
	
	}

		hedit(images=rc_name,fields='REFSPEC1',
		value=r_arc_out,add+,del-,ver-,show-,update+)
		hedit(images=bc_name,fields='REFSPEC1',
		value=b_arc_out,add+,del-,ver-,show-,update+)
		dispcor(inp=rc_name,outp=red//'_w',lin+,log-,flux+,
		same-,glob-,ignor-,confirm+,listo-,verb+)
		dispcor(inp=bc_name,outp=blue//'_w',lin+,log-,flux+,
		same-,glob-,ignor-,confirm+,listo-,verb+)

	}
}

#  Wavelength Calibration of Flux Standards	

flux_c_ld=fluxclist
flux_arcs_ld=flux_arcs

while (fscan(flux_c_ld,red,blue,cal_name)!=EOF && fscan(flux_arcs_ld,r_arc_name,b_arc_name,r_ref,b_ref)!=EOF){	
	
	flux_flag=0
	if (access(red//'_w.fits')){
		print((red),' has already been wavelength calibrated, skipping to next standard')
		flux_flag=1
		
	 }	

	if (flux_flag==0){
	total_no=obj_no+ sm_no +sm_f_no + flux_no
	b=0
	for (i=1 ; i<=total_no ; i+=1){
		if (access(r_arc_name//"_t"//b//'.fits')){
			b=b+1
		}
	}
	
	r_arc_out=r_arc_name//"_t"//b
	b_arc_out=b_arc_name//"_t"//b
	rc_name=red//'_o'
	bc_name=blue//'_o'


	
	print("Wavelength Calibrating Flux Standard ",(red)," and ",(blue))

	apall(r_arc_name,out=(r_arc_out)//'.fits',ref=(red)//'_C',recen-,trace-,
	back-,intera-,t_fun='spline3')
	
	apall(b_arc_name,out=(b_arc_out)//'.fits',ref=(blue)//'_C',recen-,trace-,
	back-,intera-,t_fun='spline3')
		

	if ((r_arc_out) == (r_ref//'_t0')){
		print ('identifying ',(r_arc_name)) 
		identify(images=r_arc_out,secti='middle line',datab='database',
		coordli='linelists$idhenear.dat',nsum=10,match=10.,
		maxfeat=75,zwidth=100.,ftype='emission',fwid=2.,cradius=4.,
		thres=20.,mins=2.,functi='legendre',order=2,niter=0,low_rej=3.,
		high_rej=3.,grow=0.,autow-,graph='stdgraph',autowri+)

		print ('identifying ',(b_arc_name)) 
		identify(images=b_arc_out,secti='middle line',datab='database',
		coordli='linelists$idhenear.dat',nsum=10,match=10.,
		maxfeat=75,zwidth=100.,ftype='emission',fwid=2.,cradius=4.,
		thres=20.,mins=2.,functi='legendre',order=2,niter=0,low_rej=3.,
		high_rej=3.,grow=0.,autow-,graph='stdgraph',autowri+)
		

	}

	if ((r_arc_out) != (r_ref//'_t0')){
	print ('reidentifying ',(r_arc_name), ' using reference ',(r_ref))
 
		reidentify(images=r_arc_out,refer=r_ref//'_t0',interac+,
		secti='middle line',newaps+,overr+,refit+,trace-,step=10,
		nsum=10,shift=0.,nlost=0,crad=5.,thresh=10.,add-,
		coord='linelists$idhenear.dat',match=10.,maxf=50,mins=2.,
		datab='database',logf='logfile',verb-,graph='stdgraph',ans+)

	print ('reidentifying ',(b_arc_name), ' using reference ',(b_ref))
 
		reidentify(images=b_arc_out,refer=b_ref//'_t0',interac+,
		secti='middle line',newaps+,overr+,refit+,trace-,step=10,
		nsum=10,shift=0.,nlost=0,crad=5.,thresh=10.,add-,
		coord='linelists$idhenear.dat',match=10.,maxf=50,mins=2.,
		datab='database',logf='logfile',verb-,graph='stdgraph',ans+)
	
	}

		hedit(images=rc_name,fields='REFSPEC1',
		value=r_arc_out,add+,del-,ver-,show-,update+)
		hedit(images=bc_name,fields='REFSPEC1',
		value=b_arc_out,add+,del-,ver-,show-,update+)
		dispcor(inp=rc_name,outp=red//'_w',lin+,log-,flux+,
		same-,glob-,ignor-,confirm+,listo-,verb+)
		dispcor(inp=bc_name,outp=blue//'_w',lin+,log-,flux+,
		same-,glob-,ignor-,confirm+,listo-,verb+)

	}	
}

sm_f_std_ld='sm_c_f.list'
flux_c_ld=fluxclist

# Sky Absorption Correction for Flux Standards

print ('Sky Absorption Correction for Flux Standards')

while (fscan(sm_f_std_ld,r_sm,b_sm)!=EOF && fscan(flux_c_ld,r_flux,b_flux,cal_name)!=EOF){
		if (access(r_sm//'_c.fits')){
		if (access(b_sm//'_c.fits')){
		print ((r_sm),' and ',(b_sm),' have already had their continuum fitted')
		goto end_cont_fit
		}
	}
	
	print ('Fitting the continuum to smooth spectrum standard ',(r_sm))

	continuum(input=r_sm//'_w',output=r_sm//'_c',type='ratio',replace-,
	wave+,logsc-,overr-,liston-,logfile='logfile',inter+,sample='*',
	naver=1,funct='spline3',order=3,low_rej=2.5,high_rej=3.0,niter=10,
	grow=1.,markre+,graphi='stdgraph')

	print ('Fitting the continuum to smooth spectrum standard ',(b_sm))

	continuum(input=b_sm//'_w',output=b_sm//'_c',type='ratio',replace-,
	wave+,logsc-,overr-,liston-,logfile='logfile',inter+,sample='*',
	naver=1,funct='spline3',order=3,low_rej=2.5,high_rej=3.0,niter=10,
	grow=1.,markre+,graphi='stdgraph')
	
end_cont_fit:

	print ('Correcting smooth spectrum template ',(r_sm),' for Halpha absorption and noise')


		if (access(r_sm//'_ca.fits')){
		if (access(b_sm//'_ca.fits')){
		print ((r_sm),' and ',(b_sm),' have already been made into templates')
		goto end_temp
		}
	}
r_template1:	
	splot(images=r_sm//'_c',line=1,band=1,units='',
	options='auto wreset xydraw hist',xmin=INDEF,xmax=INDEF,ymin=INDEF,
	ymax=INDEF,save_fi='splot.log',graphic='stdgraph',functi='spline3',
	order=11,low_rej=3,high_rej=3,niter=10,grow=1,markrej+,
	new_image=r_sm//'_ca',overwr-)
	
	if (!access(r_sm//'_ca.fits')){
		print ('output file not produced, repeating previous step - please type i when finished correcting template')
		goto r_template1}
	print ('Correcting smooth spectrum template ',(b_sm),' for Hbeta absorption and noise')
	
b_template1:
	splot(images=b_sm//'_c',line=1,band=1,units='',
	options='auto wreset xydraw hist',xmin=INDEF,xmax=INDEF,ymin=INDEF,
	ymax=INDEF,save_fi='splot.log',graphic='stdgraph',functi='spline3',
	order=11,low_rej=3,high_rej=3,niter=10,grow=1,markrej+,
	new_image=b_sm//'_ca',overwr-)
	
	if (!access(b_sm//'_ca.fits')){
		print ('output file not produced, repeating previous step - please type i when finished correcting template')
		goto b_template1}

end_temp:
	
	if (access(r_flux//'_d.fits')){
	if (access(b_flux//'_d.fits')){
	print ((r_flux),' and ',(b_flux),' have already been corrected for sky absorption')
	goto end_abs_corr
	}
	}
	
	print('Correcting flux standard ',(r_flux),' for sky absorption') 
		
	sarith(input1=r_flux//'_w',op='/',input2=r_sm//'_ca',
	outp=r_flux//'_d',w1=INDEF,w2=INDEF,apert='',bands='',
	clobber-,ignoreap+,reverse-,offset=0,renum-,merge-,rebin+,
	errval=0.,verb-)
	
	print('Correcting flux standard ',(b_flux),' for sky absorption') 
		
	sarith(input1=b_flux//'_w',op='/',input2=b_sm//'_ca',
	outp=b_flux//'_d',w1=INDEF,w2=INDEF,apert='',bands='',
	clobber-,ignoreap+,reverse-,offset=0,renum-,merge-,rebin+,
	errval=0.,verb-)
	
end_abs_corr:
}

# Making file of flux and smooth spectrum standards

flux_c_ld=fluxclist
sm_c_ld='sm_c.list'
if (access("stand.list")){
		delete 'stand.list'
}
if(access('r_sm_temp.list')){
		delete 'r_sm_temp.list'
}
if(access('b_sm_temp.list')){
		delete 'b_sm_temp.list'
}




while (fscan(flux_c_ld,r_flux,b_flux,cal_name)!=EOF && fscan(sm_c_ld,r_sm,b_sm)!=EOF){
	print((r_flux),' ',(b_flux),' ',(r_sm),' ',(b_sm), ' ',(cal_name), >> 'stand.list')
	print((r_sm), >> 'r_sm_temp.list')
	print((b_sm), >> 'b_sm_temp.list')
}


# Correcting for Halpha and Hbeta absorption
# Halpha - 6563 Angstroms
# Hbeta - 4861 Angstroms



# Flux calibration and Sky Absorption Correction

print ('Sky Absorption Correction for Objects')

stand_ld ='stand.list'
images_c_ld='image_c.list'
r_sm=''
b_sm=''

while (fscan(stand_ld,r_flux,b_flux,r_sm,b_sm,cal_name)!=EOF && fscan(images_c_ld,red,blue)!=EOF){
	
	print ('Fitting the continuum to smooth spectrum standard',(r_sm))
	
	if (access(r_sm//'_c.fits')){
		if (access(b_sm//'_c.fits')){
		print ((r_sm),' and ',(b_sm),' have already had their continuum fitted')
		goto end_cont_fit2
		}
	}
	continuum(input=r_sm//'_w',output=r_sm//'_c',type='ratio',replace-,
	wave+,logsc-,overr-,liston-,logfile='logfile',inter+,sample='*',
	naver=1,funct='spline3',order=3,low_rej=2.5,high_rej=3.0,niter=10,
	grow=1.,markre+,graphi='stdgraph')

	print ('Fitting the continuum to smooth spectrum standard',(b_sm))

	continuum(input=b_sm//'_w',output=b_sm//'_c',type='ratio',replace-,
	wave+,logsc-,overr-,liston-,logfile='logfile',inter+,sample='*',
	naver=1,funct='spline3',order=3,low_rej=2.5,high_rej=3.0,niter=10,
	grow=1.,markre+,graphi='stdgraph')

end_cont_fit2:	
	
	print ('Correcting smooth spectrum template ',(r_sm),' for Halpha absorption and noise')
	if (access(r_sm//'_ca.fits')){
		if (access(b_sm//'_ca.fits')){
		print ((r_sm),' and ',(b_sm),' have already been made into templates')
		goto end_temp2
		}
	}
r_template:	
	splot(images=r_sm//'_c',line=1,band=1,units='',
	options='auto wreset xydraw hist',xmin=INDEF,xmax=INDEF,ymin=INDEF,
	ymax=INDEF,save_fi='splot.log',graphic='stdgraph',functi='spline3',
	order=11,low_rej=3,high_rej=3,niter=10,grow=1,markrej+,
	new_imag=r_sm//'_ca',overwr-)
	
	if (!access(r_sm//'_ca.fits')){
		print ('output file not produced, repeating previous step - please type i when finished correcting template')
		goto r_template}
	print ('Correcting smooth spectrum template ',(b_sm),' for Hbeta absorption and noise')
	
b_template:
	splot(images=b_sm//'_c',line=1,band=1,units='',
	options='auto wreset xydraw hist',xmin=INDEF,xmax=INDEF,ymin=INDEF,
	ymax=INDEF,save_fi='splot.log',graphic='stdgraph',functi='spline3',
	order=11,low_rej=3,high_rej=3,niter=10,grow=1,markrej+,
	new_imag=b_sm//'_ca',overwr-)
	
	if (!access(b_sm//'_ca.fits')){
		print ('output file not produced, repeating previous step - please type i when finished correcting template')
		goto b_template}


	
end_temp2:

	for (i=1;i<=wid_no;i+=1){
		print('Correcting image ',(red),' aperture ',(i),' for sky absorption') 
		
		if (access(red//'_d'//i//'.fits')){
		if (access(blue//'_d'//i//'.fits')){
	print ((red),' and ',(blue),' have already been corrected for sky absorption')
	goto end_abs_corr2
	}
	}

		sarith(input1=red//'_w'//i,op='/',input2=r_sm//'_ca',
		outp=red//'_d'//i,w1=INDEF,w2=INDEF,apert='',bands='',
		clobber-,ignoreap+,reverse-,offset=0,renum-,merge-,rebin+,
		errval=0.,verb-)

	
		print('Correcting image ',(blue),' aperture ',(i),' for sky absorption') 
		
		sarith(input1=blue//'_w'//i,op='/',input2=b_sm//'_ca',
		outp=blue//'_d'//i,w1=INDEF,w2=INDEF,apert='',bands='',
		clobber-,ignoreap+,reverse-,offset=0,renum-,merge-,rebin+,
		errval=0.,verb-)
end_abs_corr2:
	}

	print('Flux calibrating object ',(red))


	if (access(r_flux//'_std.fits')){
		goto flux_stage1
		}

	if (access(r_flux//'_da.fits')){
		goto flux_stage0
		}
flux_loop1:

	print('correcting flux calibrator for Halpha absorption, type x to extrapolate over absorption, and i to save')

	splot(images=r_flux//'_d',line=1,band=1,units='',
	options='auto wreset xydraw hist',xmin=INDEF,xmax=INDEF,ymin=INDEF,
	ymax=INDEF,save_fi='splot.log',graphic='stdgraph',functi='spline3',
	order=11,low_rej=3,high_rej=3,niter=10,grow=1,markrej+,
	new_imag=r_flux//'_da',overwr-)

flux_stage0:
	
	if (!access(r_flux//'_da.fits')){
		print('please press i to save image')
		goto flux_loop1
		}
		
	imgets(r_flux//'_d',param='exposed')
	r_exposed=real(imgets.value)
	imgets(red//'_w1',param='exposed')
	r_obj_exposed=real(imgets.value)

	standard(input=r_flux//'_da',outp=r_flux//'_std',sames+,beam-,
	aperture='',bandwid=INDEF,bandsep=INDEF,fnuzero=3.680E-20,
	extinct='onedstds$ssoextinct.dat',caldir=(cal_dir),
	star_nam=(cal_name),inter+,
	exptime=(r_exposed),observa=')_.observatory',graphic='stdgraph')
	
flux_stage1:

	if (access(r_flux//'_sens.0001.fits')){
		goto flux_stage2
	}
	
	sensfunc(standard=r_flux//'_std',sens=r_flux//'_sens',ignor-,
	newextinc='extinct.dat',observa=')_.observatory',functio='spline3',
	order=6,inter+)

flux_stage2:

	for (i=1;i<=wid_no;i+=1){
		if (access(red//'_f'//i//'.fits')){
		goto flux_stage3}

		calibrate(input=red//'_d'//i,outp=red//'_f'//i,
		sens=r_flux//'_sens.0001',exptime=(r_obj_exposed),
		extinction='onedstds$ssoextinct.dat',
		flux+,extinct+,observa=')_.observatory',ignorea+,fnu-)
flux_stage3:					
	}

	print('Flux calibrating object ',(blue))

	
	if (access(b_flux//'_da.fits')){
		goto flux_stage3_5
		}
flux_loop2:

	print('correcting flux calibrator for Hbeta absorption, type x to extrapolate over absorption, and i to save')

	splot(images=b_flux//'_d',line=1,band=1,units='',
	options='auto wreset xydraw hist',xmin=INDEF,xmax=INDEF,ymin=INDEF,
	ymax=INDEF,save_fi='splot.log',graphic='stdgraph',functi='spline3',
	order=11,low_rej=3,high_rej=3,niter=10,grow=1,markrej+,
	new_imag=b_flux//'_da',overwr-)

flux_stage3_5:

	if (!access(r_flux//'_da.fits')){
		print('please press i to save image')
		goto flux_loop2
		}

	imgets(b_flux//'_d',param='exposed')
	b_exposed=real(imgets.value)
	imgets(blue//'_w1',param='exposed')
	b_obj_exposed=real(imgets.value)

	if (access(b_flux//'_std.fits')){
		goto flux_stage4}

	standard(input=b_flux//'_da',outp=b_flux//'_std',sames+,beam-,
	aperture='',bandwid=INDEF,bandsep=INDEF,fnuzero=3.680E-20,
	extinct='onedstds$ssoextinct.dat',caldir=(cal_dir),inter+,
	star_nam=(cal_name),inter+,
	exptime=(b_exposed),observa=')_.observatory',graphic='stdgraph')
	
flux_stage4:
	
	if (access(b_flux//'_sens.0001.fits')){
		goto flux_stage5}
	
	sensfunc(standard=b_flux//'_std',sens=b_flux//'_sens',ignor-,
	newextinc='extinct.dat',observa=')_.observatory',functio='spline3',
	order=6,inter+)

flux_stage5:

	for (i=1;i<=wid_no;i+=1){

		if (access(blue//'_f'//i//'.fits')){
		goto flux_stage6 }

		calibrate(input=blue//'_d'//i,outp=blue//'_f'//i,
		sens=b_flux//'_sens.0001',exptime=(b_obj_exposed),
		extinction='onedstds$ssoextinct.dat',
		flux+,extinct+,observa=')_.observatory',ignorea+,fnu-)
					
flux_stage6:	
	}
	
}

print ('Finished!')
print ('WooHoo!')

cleanup:
if (access('b_flux.list')){
	delete 'b_flux.list'}
if (access('r_flux.list')){
	delete 'r_flux.list'}
if (access('b_sm.list')){
	delete 'b_sm.list'}
if (access('r_sm.list')){
	delete 'r_sm.list'}
if (access('b_sm_temp.list')){
	delete 'b_sm_temp.list'}
if (access('barc.list')){
	delete 'barc.list'}
if (access('rarc.list')){
	delete 'rarc.list'}
if (access('bflux.list')){
	delete 'bflux.list'}
if (access('rflux.list')){
	delete 'rflux.list'}
if (access('bflux_c.list')){
	delete 'bflux_c.list'}
if (access('rflux_c.list')){
	delete 'rflux_c.list'}
if (access('bimage_c.list')){
	delete 'bimage_c.list'}
if (access('rimage_c.list')){
	delete 'rimage_c.list'}
if (access('blue_ctemp.list')){
	delete 'blue_ctemp.list'}
if (access('red_ctemp.list')){
	delete 'red_ctemp.list'}
if (access('bsm.list')){
	delete 'bsm.list'}
if (access('rsm.list')){
	delete 'rsm.list'}
if (access('cal_name.list')){
	delete 'cal_name.list'}
if (access('fbcomb.list')){
	delete 'fbcomb.list'}
if (access('frcomb.list')){
	delete 'frcomb.list'}
if (access('flux_c.list')){
	delete 'flux_c.list'}
if (access('image_c.list')){
	delete 'image_c.list'}
if (access('sm_bcomb.list')){
	delete 'sm_bcomb.list'}
if (access('sm_rcomb.list')){
	delete 'sm_rcomb.list'}
if (access('sm_c.list')){
	delete 'sm_c.list'}
if (access('stand.list')){
	delete 'stand.list'}
if (access('std.list')){
	delete 'std.list'}
if (access('rsm_t.list')){
		delete 'rsm_t.list'}
if (access('bsm_t.list')){
		delete 'bsm_t.list'}
if (access('rsm_f.list')){
		delete 'rsm_f.list'}
if (access('bsm_f.list')){
		delete 'bsm_f.list'}

end