#pragma TextEncoding = "UTF-8"
#pragma version=5.10
#pragma rtGlobals=3
#pragma IgorVersion=9

//The procedure was developed by Dr. Gregor Germer within a research project of Prof. Dr. Eckart Rühl (Physical Chemistry, Freie Universität Berlin, Berlin, Germany) funded by the German Research Foundation (DFG), grant no: RU 420/12-1.
//
//It is required to cite the following reference when this procedure is used for scientific work and publications: G. Germer, T. Ohigashi, H. Yuzawa, N. Kosugi, R. Flesch, F. Rancan, A. Vogt, and E. Rühl, “Improved Skin Permeability after Topical Treatment with Serine Protease: Probing the Penetration of Rapamycin by Scanning Transmission X-Ray Microscopy”, ACS Omega 6, 12213-12222 (2021). doi.org/10.1021/acsomega.1c01058. 
//
//Many thanks for the help from the wavemetrics forum, especially to the users tony, chozo, jjweimer, oleytken, KZarzana, KurtB, johnweeks (WM) and aclight (WM). Without their help this extended version wouldn’t achieve the current state.
//Thanks to my students for beta testing
//
//changelog:
//ggermer - 5.10: add support for 4-dim. wave-to-fit waves - Please be careful, I don't have a measured dataset to check, the check on a generated dataset seems to be ok.
//ggermer - 5.07: saves the note of the to-be-fitted-wave in the information wave, shorter names for LCM results & name of to-be-fitted-wave saved as note in LCM_coef & LCM_sigma
//ggermer - 5.06: LCM panel position is relative to the middle of the active screen / the screen where the mouse is (Macintosh) or Igors windows (Windows), Waves panel relative to LCM panel, fixedSize=1 for the LCMpanel, Panel position optimation using Windows
//ggermer - 5.05: restart bug fixes, correction of start/end for x-shift pay attention to the initial value
//ggermer - 5.04: MT optimization 1-dim. waves, uses version numbering compatible with the pragma version, mask_red restart bug fixed, change of units, bug fix 1-dim. fit_wave, waves w_coef & w_sigma & w_references now killed
//ggermer - 5.03: save LCM configration files, bug fix for wave with names containing space (e.g. MPF results)
//ggermer - 5.02: single energy fit for all reference waves together, restart LCM using information wave, restart LCM using old mask wave, deactivation of point-by-point-fits, bug fix restart LCM after aborted LCM (MT problem), generate fit wave after MT process
//ggermer - 5.01: bug fix for hold function, bug fix LCM panel (green border), using of double float waves, bug fix for 3-dim. waves: (x,1,z) and (1,y,z)
//ggermer - 5.00: completely rewritten version, now without the use quick and dirty methods, new layout for panel, drop-down menu for panel, help text for panel, new spectrums panel
//ggermer - 4.xx: multiple energy shift for all reference waves, all-at-ones funcfit, results of 1 dim. waves saved as global variables, buttons for fast changing of initial guess & min & max, LCM funcfit uses /NWOK, bug fix all-at-ones funcfit
//ggermer - 3.xx: average spectrum in panel, bug fix open LCM panel if panel already exists, save LCM results in a subfolder of the input, extract LCM results in multiple 2-dim. waves
//ggermer - 0-2:  new panel, mask regions as red region(s) in average spectrum, multiple bug fixes

Menu "Analysis"
	submenu "Packages"
		submenu "Linear-combination modeling"
			"Start LCM (new dataset).", f_LCM_start()
			"Restart LCM (dataset of information wave from Data Browser, no mask wave).", f_LCM_restart(0)
			"Restart LCM (dataset of information wave from Data Browser, using mask wave).", f_LCM_restart(1)
			"Calculate relative sigma of LCM sigma wave from Data Browser.", f_sigma_to_sigma_rel()
			"Release ThreadGroup for aborted LCM using Igor Pro's beach ball.", f_ThreadGroupRelease_all()
		end
	end
end

//----------------------------------------------------------------------------

function f_LCM_beep()
	variable v_beep=0 //0 = no beep, 1 = beep
	return v_beep
end


function/S f_LCM_version()
	string s_version=num2str(ProcedureVersion("f_LCM_start"))
	return s_version
end


function f_LCM_start()
	if (WinType("LCMpanel"))	//does an active LCM exists?
		DoWindow/F LCMpanel
		return 0
	else
	//export folder generation
		f_ThreadGroupRelease_all()
		setDataFolder root:
		if (datafolderExists("root:Packages")==0)
			newdataFolder root:Packages
		endif		
		if (datafolderExists("root:Packages:LCM")==0)
			newdataFolder root:Packages:LCM
		else		//remove old files
		 	killdataFolder root:Packages:LCM
		 	newdatafolder root:Packages:LCM
		endif
		setDataFolder root:Packages:LCM										
		string testfolder = "root:Packages:LCM"
	//transfer waves
		setdataFolder $testfolder
		make/D/O/WAVE/N=1 reference_waves
		make/D/O/WAVE/N=1 reference_waves_fit
		make/D/O/WAVE/N=1 data_wave
		make/D/O/N=10 lcm_variables //[0] = v_lcm_type: 1 = all-at-ones 2 = point-by-point , [1] = x_shift_done: 1 = no 2 = yes , [2] = quantity of waves, [3] = v_LCM_offset, [4] = v_LCM_delta, [5] = v_LCM_max , [6] = fitTol, [7] = start, [8] = end, [9] = dimensions of sample_wave
	//import of sample wave
		setdataFolder root:
		variable index,index2
		variable v_num_of_waves=0
		string s_nameofwave
		string s_wavepath
		string s_promptStr = "Choose the sample wave (1, 2 or 3 dimensions)"
		CreateBrowser prompt=s_promptStr, showWaves=1, showVars=0, showStrs=0
		if (V_flag==0)
			print "LCM: Abort by user."
			killdataFolder/Z root:Packages:LCM
			return 0
		endif
		s_wavepath=removeending(s_BrowserList,";")
		if (waveexists($s_wavepath)==0)
			print "LCM: No sample wave."
			killdataFolder/Z root:Packages:LCM
			return 0
		endif	
		wave w_sample=$s_wavepath
		variable v_rows=dimsize(w_sample,0)
		variable v_columns=dimsize(w_sample,1)
		variable v_layer=dimsize(w_sample,2)
		variable v_chunks=dimsize(w_sample,3)
		string s_exportname
		variable v_transformed_2_3D_loc=0
		variable v_sampleoffset, v_sampledelta, v_rows_references
		variable v_layer_sample
		setdataFolder $testfolder
		if (v_chunks>1)
				v_sampleoffset = dimoffset(w_sample,3)
				v_sampledelta = dimdelta(w_sample,3)
				make/D/N=(v_chunks)/O energyscale=v_sampleoffset+x*v_sampledelta
				v_rows_references=v_chunks
				lcm_variables[9]=4
				v_layer_sample=dimsize(w_sample,3)		
		else
			if (v_layer>0)
				v_sampleoffset = dimoffset(w_sample,2)
				v_sampledelta = dimdelta(w_sample,2)
				make/D/N=(v_layer)/O energyscale=v_sampleoffset+x*v_sampledelta
				v_rows_references=v_layer
				lcm_variables[9]=3
				v_layer_sample=dimsize(w_sample,2)
			else
				if (v_columns>0)
					v_sampleoffset = dimoffset(w_sample,1)
					v_sampledelta = dimdelta(w_sample,1)
					make/D/N=(v_columns)/O energyscale=v_sampleoffset+x*v_sampledelta
					v_rows_references=v_columns
					lcm_variables[9]=2
					v_layer_sample=dimsize(w_sample,1)
				else
					v_sampleoffset = dimoffset(w_sample,0)
					v_sampledelta = dimdelta(w_sample,0)	
					make/D/N=(v_rows)/O energyscale=v_sampleoffset+x*v_sampledelta	
					v_rows_references=v_rows	
					lcm_variables[9]=1
					v_layer_sample=dimsize(w_sample,0)
				endif
			endif
		endif
		data_wave[0]=$s_wavepath		
	//import of reference waves
		s_promptStr = "Choose all references (1 dimension) wave for linear combination."
		setdataFolder root:
		CreateBrowser executeMode=2, prompt=s_promptStr, showWaves=1, showVars=0, showStrs=0
		if (V_flag==0)
			print "LCM: Abort by user."
			killdataFolder/Z root:Packages:LCM
			return 0
		endif		
		setdataFolder $testfolder
		variable index_max=itemsinList(S_BrowserList,";")
		variable v_position
		variable wave_index=0
		variable v_dim_wave
		string s_test
		string s_order=""
		for (index=0;index<index_max;index+=1)
			s_test=StringFromList(index, S_BrowserList, ";" )
			if (waveexists($s_test)==1)
				if (wavetype($s_test,1)==1)
					if (wave_index>0)
						insertPoints v_dim_wave,1,reference_waves
						insertPoints v_dim_wave,1,reference_waves_fit
					endif
					s_order+=s_test+";"
					s_nameofwave=testfolder+":reference"+num2str(wave_index)
					make/D/O w_inputwave
					killwaves w_inputwave
					wave w_inputwave=$s_test
					make/D/N=(v_rows_references)/O w_scaledwave
					for (index2=0;index2<v_rows_references;index2+=1)				//get the same scale
						v_position=energyscale[index2]
						v_position=((v_position - DimOffset(w_inputwave,0)) / DimDelta(w_inputwave,0))
						w_scaledwave[index2]=w_inputwave[v_position]
					endfor
					setscale/P x,v_sampleoffset,v_sampledelta,w_scaledwave
					duplicate w_scaledwave $s_nameofwave
					killwaves w_scaledwave
					reference_waves[wave_index]=$s_test
					reference_waves_fit[wave_index]=$s_nameofwave
					v_dim_wave=dimsize(reference_waves,0)
					wave_index+=1
				endif
			endif
		endfor
		lcm_variables[2]=wave_index
		v_num_of_waves=wave_index
		setdataFolder $testfolder
		if (wave_index<1)
			print "LCM: No reference waves."
			return 0	
		endif
	
	//generate MT fit wave
		make/D/O/N=(v_num_of_waves,v_layer_sample) w_references_MT
		for (index=0;index<v_num_of_waves;index+=1)
			s_nameofwave=testfolder+":reference"+num2str(index)
			if (waveexists($s_nameofwave)==0)
				print "Error 4"
				print s_nameofwave
				return 0
			endif
			wave w_inputwave=$s_nameofwave
			for (index2=0;index2<v_layer_sample;index2+=1)
				w_references_MT[index][index2]=w_inputwave[index2]
			endfor
		endfor
		setscale/P y,v_sampleoffset,v_sampledelta,w_references_MT	
		
		
	//generate preview wave
		//print v_layer,v_columns,v_rows
		variable r_index, c_index, l_index
		if (v_chunks>1)
			make/N=(v_rows,(v_columns*v_layer),v_chunks)/FREE w_3dim
			variable ch_index
			for (r_index=0;r_index<v_rows;r_index++)
				for (c_index=0;c_index<v_columns;c_index++)
					for (l_index=0;l_index<v_layer;l_index++)
						for (ch_index=0;ch_index<v_chunks;ch_index++)
							w_3dim[r_index][(c_index*v_layer)+l_index][ch_index]=w_sample[r_index][c_index][l_index][ch_index]
						endfor
					endfor
				endfor	
			endfor
			make/D/O W_ISBeamAvg
			make/D/O W_ISBeamMax
			make/D/O W_ISBeamMin
			ImageStats/BEAM/Q w_3dim 
			setscale/P x,dimoffset(w_sample,3),dimdelta(w_sample,3),W_ISBeamAvg
			duplicate/O W_ISBeamAvg root:Packages:LCM:w_example
			killwaves W_ISBeamAvg
			killwaves W_ISBeamMax
			killwaves W_ISBeamMin			
		else	
			if (v_layer>0)
				if (v_columns>1&&v_rows>1)
					make/D/O W_ISBeamAvg
					make/D/O W_ISBeamMax
					make/D/O W_ISBeamMin
					ImageStats/BEAM/Q w_sample 
					setscale/P x,dimoffset(w_sample,2),dimdelta(w_sample,2),W_ISBeamAvg
					make/D/O root:Packages:LCM:w_example
					killwaves root:Packages:LCM:w_example
					duplicate W_ISBeamAvg root:Packages:LCM:w_example
					killwaves W_ISBeamAvg
					killwaves W_ISBeamMax
					killwaves W_ISBeamMin	
				else
					if (v_columns==1&&v_rows==1)
						make/D/O/N=(v_layer) w_build
						matrixOP/O w_build=beam(w_sample,0,0)
						setscale/P x,dimoffset(w_sample,2),dimdelta(w_sample,2),w_build
						duplicate/O w_build root:Packages:LCM:w_example								
					else
						if (v_rows>1)
							make/D/O/N=(v_layer) w_build
							matrixOP/O w_build=sumCols(w_sample)
							make/D/O/N=(v_layer) w_build2
							matrixOP/O w_build2=beam(w_build,0,0)
							setscale/P x,dimoffset(w_sample,2),dimdelta(w_sample,2),w_build2
							duplicate/O w_build2 root:Packages:LCM:w_example					
						else
							make/D/O/N=(v_layer) w_build
							matrixOP/O w_build=sumRows(w_sample)
							make/D/O/N=(v_layer) w_build2
							matrixOP/O w_build2=beam(w_build,0,0)
							setscale/P x,dimoffset(w_sample,2),dimdelta(w_sample,2),w_build2
							duplicate/O w_build2 root:Packages:LCM:w_example									
						endif
				
					endif

				endif		
			else
				if (v_columns>0)
					make/D/O/N=(v_columns) w_build
					matrixOP/O w_build=sumCols(w_sample)
					make/D/O/N=(v_columns) w_build2
					for (index=0;index<v_columns;index++)
						w_build2[index]=w_build[0][index]
					endfor
					setscale/P x,dimoffset(w_sample,1),dimdelta(w_sample,1),w_build2
					duplicate/O w_build2 root:Packages:LCM:w_example	
				else
					duplicate/O w_sample root:Packages:LCM:w_example			
				endif
			endif
		endif
	//generate panel
		//lcm_variables //[0] = v_lcm_type: 1 = all-at-ones 2 = point-by-point , [1] = x_shift_done: 1 = no 2 = yes , [2] = quantity of waves, [3] = v_LCM_offset, [4] = v_LCM_delta, [5] = v_LCM_max, [6] = fitTol, [7] = start, [8] = end, [9] = dimensions of sample_wave
		make/D/O/N=(v_rows_references) root:Packages:LCM:w_mask=1
		make/D/O/N=(v_rows_references) root:Packages:LCM:w_mask_red=NaN
		wave w_mask=root:Packages:LCM:w_mask
		wave w_mask_red=root:Packages:LCM:w_mask_red		
		setscale/P x,v_sampleoffset,v_sampledelta,w_mask 
		setscale/P x,v_sampleoffset,v_sampledelta,w_mask_red
		variable v_x_shift_abs_max=floor(((v_sampleoffset+v_sampledelta*(v_rows_references-1))-v_sampleoffset)*0.45)
		variable v_x_shift_start_max=floor(((v_sampleoffset+v_sampledelta*(v_rows_references-1))-v_sampleoffset)*0.2)
		lcm_variables[3]=v_sampleoffset
		lcm_variables[4]=v_sampledelta
		lcm_variables[5]=(v_sampleoffset+v_sampledelta*(v_rows_references-1))
	//screen resolution	
		string s_igor_info=IgorInfo(0)
		variable v_info_elements=itemsinList(s_igor_info,";")
		string s_screen_info=replaceString("RECT=",stringfromList(v_info_elements-1,s_igor_info,";"),"")
		variable v_screen1_width=abs(str2num(stringfromList(3,s_screen_info,","))-str2num(stringfromList(1,s_screen_info,","))) //relative position
		variable v_screen1_length=abs(str2num(stringfromList(4,s_screen_info,","))-str2num(stringfromList(2,s_screen_info,",")))
		variable v_sc1_width_mp=round(v_screen1_width/2)
		variable v_sc1_length_mp=round(v_screen1_length/2)
		if (stringmatch(IgorInfo(2),"Windows")==1)
			GetWindow kwFrameInner, wsize
			v_sc1_width_mp=round((V_right-V_left)/2) //relative position
			v_sc1_length_mp=round((V_bottom-V_top)/2)
		endif	
		if (stringmatch(IgorInfo(2),"Macintosh")==1)
			GetMouse
			string s_Igorinfo = IgorInfo(0)
			string s_str
			variable v_n_screens=NumberByKey("NSCREENS", s_Igorinfo)
			//print v_n_screens
			variable v_i, v_screen_left, v_screen_top, v_screen_right, v_screen_bottom
			make/FREE/O/N=(v_n_screens,4) w_monitors
			for (v_i=1;v_i<=v_n_screens;v_i++)
				s_str = StringByKey("SCREEN"+num2str(v_i), s_Igorinfo)
				sscanf s_str, "DEPTH=%*d,RECT=%d,%d,%d,%d", v_screen_left, v_screen_top, v_screen_right, v_screen_bottom
				//print v_screen_left, v_screen_top, v_screen_right, v_screen_bottom
				w_monitors[v_i-1][0]=v_screen_left
				w_monitors[v_i-1][1]=v_screen_right
				w_monitors[v_i-1][2]=v_screen_top
				w_monitors[v_i-1][3]=v_screen_bottom
			endfor
			variable v_active_screen=1
			for (v_i=1;v_i<=v_n_screens;v_i++)
				if ((v_left>=w_monitors[v_i-1][0])&&(v_left<=w_monitors[v_i-1][1])&&(v_top>=w_monitors[v_i-1][2])&&(v_top<=w_monitors[v_i-1][3]))
					v_active_screen=v_i 
				endif
			endfor
			//print v_active_screen
			s_str = StringByKey("SCREEN"+num2str(v_active_screen), s_Igorinfo)
			sscanf s_str, "DEPTH=%*d,RECT=%d,%d,%d,%d", v_screen_left, v_screen_top, v_screen_right, v_screen_bottom
			//print v_screen_left, v_screen_right
			//print v_screen_top, v_screen_bottom
			v_sc1_width_mp=round((v_screen_right+v_screen_left)/2) //absolute position, negative values possible
			v_sc1_length_mp=round((v_screen_bottom+v_screen_top)/2)
		endif	
		//print v_sc1_width_mp,	v_sc1_length_mp
	//window for data input
		variable v_height=30*(wave_index)
		if (v_height < (75))
			v_height=(75)
		endif
		variable v_panel_left=v_sc1_width_mp-483 //165
		variable v_panel_rigth=v_sc1_width_mp+483 //(340+425+165+200+1)
		variable v_panel_top=v_sc1_length_mp-floor((v_height+110)/2) //140
		variable v_panel_bottom=v_sc1_length_mp+ceil((v_height+110)/2) //(v_height+110+140)
		make/O/N=4 root:Packages:LCM:w_panel_resolution
		wave w_panel_resolution=root:Packages:LCM:w_panel_resolution //left,right,top,bottom
		w_panel_resolution[0]=v_panel_left
		w_panel_resolution[1]=v_panel_rigth
		w_panel_resolution[2]=v_panel_top
		w_panel_resolution[3]=v_panel_bottom
		NewPanel/K=2/W=(v_panel_left,v_panel_top,v_panel_rigth,v_panel_bottom)/N=LCMpanel as "linear-combination modeling input" // /W=(165,140,(425+165),(30*(quantity_waves)+110+140))
		if (stringmatch(IgorInfo(2),"Macintosh")==1)
			MoveWindow v_panel_left,v_panel_top,-1,-1 //bug fix for negative positions on Macintosh
		endif
		ModifyPanel/W=LCMpanel fixedSize=1
		DrawRect 422,-2,(422+200),(v_height+110+140+2)
		DrawText 10	,25, "General settings"
		SetVariable V_FitTol	,pos={160	,10}	,value=_NUM:0.0001, limits={1e-10,0.01,0}, help={"Accuracy of the fit - between 0.01 and 1e-10."}
		DrawText 125	,25, "FitTol"
		SetVariable V_StartFit	,pos={270	,10}	,value=_NUM:v_sampleoffset, limits={v_sampleoffset,(v_sampleoffset+v_sampledelta*(v_rows_references-1)),0}, help={"Start value (energy) for fit."}
		DrawText 240	,25, "Start"
		SetVariable V_EndFit	,pos={365	,10}	,value=_NUM:(v_sampleoffset+v_sampledelta*(v_rows_references-1)), limits={v_sampleoffset,(v_sampleoffset+v_sampledelta*(v_rows_references-1)),0}, help={"End value (energy) for fit."}
		DrawText 340,25, "End"
		string s_LCM_sample_global=nameofwave(w_sample)
		string s_LCM_sample_local=ReplaceString("'",s_LCM_sample_global,"") //if wave starts with bad character
		string s_samplename
		variable v_start=v_sampleoffset
		variable v_end=(v_sampleoffset+v_sampledelta*(v_rows_references-1))
		if (strlen(s_LCM_sample_local)<38) //for long wavenames
			s_samplename="Reference waves for "+s_LCM_sample_local+""
		else
			s_samplename="Reference waves for "+s_LCM_sample_local[0,5]+" ... "+s_LCM_sample_local[(strlen(s_LCM_sample_local)-33),(strlen(s_LCM_sample_local)-1)]+""
		endif
		DrawText 10	,45, s_samplename
		PopupMenu popup0 pos={430,5},size={50.00,14.00},value="All-at-ones Fit;Point-by-point-Fit",help={"The All-at-ones Fit is faster than the Point-by-point-Fit."},disable=2
		PopupMenu popup1 pos={430,25},size={50.00,14.00},value="No x-shift analysis;Multiple x-shift analysis;Single x-shift analysis",help={"Better results for x-shift analysis if Initial Guess of waves >> 1e-6 (e. g. 0.5), ± = zero holds the initial value, Single x-shift analysis uses the initial x-shift of the first wave"}
		string s_start
		string s_min
		string s_max, s_x_start, s_x_var, s_x_hold
		string s_wavepos=""
		variable v_intemsinlist
		for (index=0;index<wave_index;index+=1)
			s_start="V_Start_"+num2str(index)
			s_min="V_Min_"+num2str(index)
			s_max="V_Max_"+num2str(index)
			s_x_start="V_xStart_"+num2str(index)
			s_x_var="V_xVar_"+num2str(index)
			s_wavepos=nameofwave(reference_waves[index])
			s_wavepos=ReplaceString("'",s_wavepos,"") //if wave starts with bad character
			DrawText 10	,((50+(30*index))+15), s_wavepos[0,12]
			SetVariable $s_start	,pos={175	,(50+(30*index))}	,value=_NUM:1e-6, limits={-inf,inf,0} , help={"Initial Guess for the fit, must be non-zero, e.g. 1e-6, of the wave "+getwavesDataFolder(reference_waves[index],2)+"."}
			DrawText 105	,((50+(30*index))+15), "Initial Guess"
			SetVariable $s_min	,pos={270	,(50+(30*index))}	,value=_NUM:-inf, limits={-inf,inf,0} , help={"Lower limit of the fit, usually -inf or zero"}
			DrawText 242	,((50+(30*index))+15), "Min"
			SetVariable $s_max	,pos={365	,(50+(30*index))}	,value=_NUM:+inf, limits={-inf,inf,0} , help={"Upper limit of the fit, usually +inf"}
			DrawText 340,((50+(30*index))+15), "Max"
			SetVariable $s_x_start	,pos={500	,(50+(30*index))}	,value=_NUM:0, limits={-v_x_shift_start_max,v_x_shift_start_max,0} , help={"Initial Guess, can be zero"}
			DrawText 430	,((50+(30*index))+15), "Initial x-shift"
			SetVariable $s_x_var	,pos={565	,(50+(30*index))}	,value=_NUM:0, limits={-v_x_shift_abs_max,v_x_shift_abs_max,0} , help={"Range for energyshift, e.g. "+num2str(10*v_sampledelta)+", Zero is interpreted as hold value"}
			DrawText 553	,((50+(30*index))+15), "±"		
		endfor
		
		DrawText 10	,((50+(30*index))+15), "absolute offset"
		SetVariable V_start_abs	,pos={175	,(50+(30*index))}	,value=_NUM:1e-6, limits={-inf,inf,0}, help={"Initial Guess of the offset for the fit, e.g. 1e-6"}
		DrawText 105	,((50+(30*index))+15), "Initial Guess"
		SetVariable V_min_abs	,pos={270	,(50+(30*index))}	,value=_NUM:-inf, limits={-inf,inf,0}, help={"Lower limit of the offset, usually -inf"}
		DrawText 242	,((50+(30*index))+15), "Min"
		SetVariable V_max_abs	,pos={365	,(50+(30*index))}	,value=_NUM:+inf, limits={-inf,inf,0}, help={"Upper limit of the offset, usually +inf"}
		DrawText 340,((50+(30*index))+15), "Max"
		Button B_Doit			,pos={253	,(v_height)+80}	,size={75,20}		,title="Do It" ,fColor=(257*30,257*144,257*255)		,proc=f_Button_doit_LCM, help={"Performing the fit can take several minutes. The fitted wave is "+getwavesDataFolder(w_sample,2)+"."}
		Button B_Cancel			,pos={338	,(v_height)+80}	,size={75,20}		,title="Cancel"		,proc=f_Button_cancel_LCM
		Button B_SpectraPrev,pos={10	,(v_height)+80}	,size={55,20}		,title="Waves", fColor=(52428,52428,52428)	,proc= f_Button_wave_preview_LCM, help={"Preview of the used reference waves."}
		Button B_SaveConfig,pos={75	,(v_height)+80}	,size={100,20}		,title="Save Config.",proc = f_Button_saveconfig_LCM, help={"Save the visible configuration in root:."}
		DrawText 430	,(v_height)+80+18), "Set initial Guess"
		Button B_initial_a,pos={530	,(v_height)+80}	,size={40,20} ,title="1e-6" ,proc=f_Button_initial_1em6_LCM
		Button B_initial_b,pos={575	,(v_height)+80}	,size={40,20} ,title="0.5" ,proc=f_Button_initial_0_5_LCM
		DrawText 430	,(v_height)+55+18), "Min"
		Button B_min_a,pos={452	,(v_height)+55}	,size={40,20} ,title="-inf" ,proc=f_Button_min_minf_LCM
		Button B_min_b,pos={497	,(v_height)+55}	,size={40,20} ,title="0" ,proc=f_Button_min_0_LCM	
		DrawText 548	,(v_height)+55+18), "Max"
		Button B_max_b,pos={575	,(v_height)+55}	,size={40,20} ,title="+inf" ,proc=f_Button_max_inf_LCM	
	//window part for mask
		//SetDrawLayer UserBack
		Button B_mask,pos={422+200+10.00,155.00},size={100.00,20.00},title="Mask region",proc=f_Button_Mask_NaN_LCM, help={"Spectral region does not include in the fit (red region)"}
		Button B_mask,fColor=(65535,0,0)
		Button B_unmask,pos={422+200+120.00,155.00},size={100.00,20.00},title="Unmask region",proc=f_Button_Mask_1_LCM, help={"Spectral region does include in the fit (blue region)"}
		Button B_unmask,fColor=(0,43690,65535)
		SetVariable V_B_start,pos={422+200+30.00,129.00},size={100.00,14.00},title="Region Start",value=_NUM:v_start, help={"For individual points, the start and end values must be set equal."}
		SetVariable V_B_start,limits={v_start,v_end,0}
		SetVariable V_B_end,pos={422+200+165.00,129.00},size={100.00,14.00},title="Region End",value=_NUM:v_end
		SetVariable V_B_end,limits={v_start,v_end,0}
		Button B_unmask_all,pos={422+200+230.00,155.00},size={100.00,20.00},title="Unmask all", proc=f_Button_Mask_reset_LCM, help={"The mask is reset so that all spectral regions are included in the fit."}
		Button B_unmask_all,fColor=(1,65535,33232)
		Display/L/W=(422+200+10,9,422+200+10+323,119)/HOST=# w_mask
		AppendToGraph/L w_mask_red
		AppendToGraph/R root:Packages:LCM:w_example			
		ModifyGraph hbFill=2
		ModifyGraph lSize=1.5
		SetAxis left 0,1
		ModifyGraph mode(w_mask)=7,lsize(w_mask)=0,plusRGB(w_mask)=(1,16019,65535,32768),negRGB(w_mask)=(1,16019,65535,32768)
		ModifyGraph mode(w_mask_red)=7,lsize(w_mask_red)=0,plusRGB(w_mask_red)=(65535,0,0,32768),negRGB(w_mask_red)=(65535,0,0,3)
		ModifyGraph hbFill(w_mask)=3	
		ModifyGraph rgb(w_example)=(0,0,0)
		ModifyGraph usePlusRGB=1
		ModifyGraph useNegRGB=1
		ModifyGraph mirror(bottom)=2
		ModifyGraph notation(left)=1
		ModifyGraph nticks(left)=0
		ModifyGraph noLabel(left)=1
		ModifyGraph nticks(right)=0
		ModifyGraph noLabel(right)=1
		ModifyGraph UIControl=(2^0+2^1+2^3+2^4+2^5+2^6+2^7+2^8+2^9+2^10+2^11)
		Label bottom ""
		Label left "Average"				
		ModifyGraph margin(right)=15,margin(top)=15,margin(left)=20
		RenameWindow #,G0
		SetActiveSubwindow ##		
	endif
end

function f_LCM_restart(variable v_use_mask)
	if (WinType("LCMpanel"))	//does an active LCM exists?
		DoWindow/F LCMpanel
		return 0
	endif
	f_ThreadGroupRelease_all()
	string s_wave=getbrowserSelection(0)
 	wave/T w_information=$s_wave
 	if (waveexists(w_information)==0)
 		print "Error: No wave selected."
 		return 0
 	endif
	string s_information_name=nameofwave(w_information)
	if ((grepString(s_information_name,"information")==0)&&(grepString(s_information_name,"LCM_configuration")==0))
 		print "Error: No information wave selected."
 		return 0	
	endif
	variable v_rows_information=dimsize(w_information,0)
	variable v_number_of_waves=0
	variable index
	for (index=1;index<v_rows_information;index++)
		if (grepString(w_information[index][1],"absolute Offset")==1)
 			index=inf	
		else
			v_number_of_waves++
		endif		
	endfor
	string s_reference_waves=""
	for (index=1;index<v_rows_information;index++)
		if (grepString(w_information[index][1],"absolute Offset")==1)
 			index=inf	
		else
			if (waveexists($w_information[index][1])==0)
 				print "Error: Reference wave "+w_information[index][1]+" does not exists."
 				return 0
 			else
 				s_reference_waves+=w_information[index][1]+";"
 			endif
		endif		
	endfor
	variable v_colums_information=dimsize(w_information,1)
	string s_mask_position=removeending(s_wave,"information")+"mask"
	//print v_number_of_waves
	//print s_reference_waves
	make/FREE/O/N=(v_number_of_waves+1) w_initial_guess, w_min, w_max, w_initial_xshift, w_var_xshift
	for (index=1;index<v_rows_information;index++)
		if (grepString(w_information[index][1],"absolute Offset")==1)
			w_initial_guess[index-1]=str2num(w_information[index][2])
			w_min[index-1]=str2num(w_information[index][3])
			w_max[index-1]=str2num(w_information[index][4]) 			
 			index=inf	
		else
			w_initial_guess[index-1]=str2num(w_information[index][2])
			w_min[index-1]=str2num(w_information[index][3])
			w_max[index-1]=str2num(w_information[index][4])
			if (v_colums_information>6)
				w_initial_xshift[index-1]=str2num(w_information[index][5])
				w_var_xshift[index-1]=str2num(w_information[index][6])
			else
				w_initial_xshift[index-1]=0
				w_var_xshift[index-1]=0			
			endif
		endif	
	endfor
	for (index=0;index<dimsize(w_initial_xshift,0);index++)
		if (numtype(w_initial_xshift[index])!=0)
			w_initial_xshift[index]=0
		endif
		if (numtype(w_var_xshift[index])!=0)
			w_var_xshift[index]=0
		endif		
	endfor
	string s_start_end_restart=w_information[5][0]
	if (stringmatch(s_start_end_restart,"FitTol")==1) //compatibility to older versions
		s_start_end_restart=w_information[4][0]
	endif
	s_start_end_restart=replaceString(" to ",s_start_end_restart,";")
	s_start_end_restart=replaceString("Scale: ",s_start_end_restart,"")
	//print s_start_end_restart
	variable v_start_restart=str2num(StringFromList(0,s_start_end_restart,";"))
	variable v_end_restart=str2num(StringFromList(1,s_start_end_restart,";"))
	string s_fittype_restart=w_information[5][0]
	variable v_allatones=2
	if (grepString(w_information[8][0],"All-at-ones Fit")==1)
		v_allatones=1
	endif
	variable v_xshift_restart = 1
	if (grepString(w_information[8][0],"Multiple x-shift")==1)
		v_xshift_restart=2
	endif
	if (grepString(w_information[8][0],"Single x-shift")==1)
		v_xshift_restart=3
	endif		
	variable v_fittol_restart=str2num(replaceString("FitTol: ", w_information[6][0], ""))
	//print v_start_restart
	//print v_end_restart
	//print v_fittol_restart
	//print v_allatones
	//print v_xshift_restart
	//export folder generation
		setDataFolder root:
		if (datafolderExists("root:Packages")==0)
			newdataFolder root:Packages
		endif		
		if (datafolderExists("root:Packages:LCM")==0)
			newdataFolder root:Packages:LCM
		else		//remove old files
		 	killdataFolder root:Packages:LCM
		 	newdatafolder root:Packages:LCM
		endif
		setDataFolder root:Packages:LCM										
		string testfolder = "root:Packages:LCM"
	//transfer waves
		setdataFolder $testfolder
		make/D/O/WAVE/N=1 reference_waves
		make/D/O/WAVE/N=1 reference_waves_fit
		make/D/O/WAVE/N=1 data_wave
		make/D/O/N=10 lcm_variables //[0] = v_lcm_type: 1 = all-at-ones 2 = point-by-point , [1] = x_shift_done: 1 = no 2 = yes , [2] = quantity of waves, [3] = v_LCM_offset, [4] = v_LCM_delta, [5] = v_LCM_max , [6] = fitTol, [7] = start, [8] = end, [9] = dimensions of sample_wave
	//import of sample wave
		setdataFolder root:
		variable index2
		variable v_num_of_waves=0
		string s_nameofwave
		string s_wavepath
		string s_promptStr = "Choose the sample wave (1, 2 or 3 dimensions)"
		CreateBrowser prompt=s_promptStr, showWaves=1, showVars=0, showStrs=0
		if (V_flag==0)
			print "LCM: Abort by user."
			killdataFolder/Z root:Packages:LCM
			return 0
		endif
		s_wavepath=removeending(s_BrowserList,";")
		if (waveexists($s_wavepath)==0)
			print "LCM: No sample wave."
			killdataFolder/Z root:Packages:LCM
			return 0
		endif	
		wave w_sample=$s_wavepath
		variable v_rows=dimsize(w_sample,0)
		variable v_columns=dimsize(w_sample,1)
		variable v_layer=dimsize(w_sample,2)
		variable v_chunks=dimsize(w_sample,3)
		string s_exportname
		variable v_transformed_2_3D_loc=0
		variable v_sampleoffset, v_sampledelta, v_rows_references
		variable v_layer_sample
		setdataFolder $testfolder
		if (v_chunks>1)
				v_sampleoffset = dimoffset(w_sample,3)
				v_sampledelta = dimdelta(w_sample,3)
				make/D/N=(v_chunks)/O energyscale=v_sampleoffset+x*v_sampledelta
				v_rows_references=v_chunks
				lcm_variables[9]=4
				v_layer_sample=dimsize(w_sample,3)		
		else
			if (v_layer>0)
				v_sampleoffset = dimoffset(w_sample,2)
				v_sampledelta = dimdelta(w_sample,2)
				make/D/N=(v_layer)/O energyscale=v_sampleoffset+x*v_sampledelta
				v_rows_references=v_layer
				lcm_variables[9]=3
				v_layer_sample=dimsize(w_sample,2)
			else
				if (v_columns>0)
					v_sampleoffset = dimoffset(w_sample,1)
					v_sampledelta = dimdelta(w_sample,1)
					make/D/N=(v_columns)/O energyscale=v_sampleoffset+x*v_sampledelta
					v_rows_references=v_columns
					lcm_variables[9]=2
					v_layer_sample=dimsize(w_sample,1)
				else
					v_sampleoffset = dimoffset(w_sample,0)
					v_sampledelta = dimdelta(w_sample,0)	
					make/D/N=(v_rows)/O energyscale=v_sampleoffset+x*v_sampledelta	
					v_rows_references=v_rows	
					lcm_variables[9]=1
					v_layer_sample=dimsize(w_sample,0)
				endif
			endif
		endif
		data_wave[0]=$s_wavepath			

		variable index_max=itemsinList(s_reference_waves,";")
		variable v_position
		variable wave_index=0
		variable v_dim_wave
		string s_test
		string s_order=""
		for (index=0;index<index_max;index+=1)
			s_test=StringFromList(index, s_reference_waves, ";" )
			if (waveexists($s_test)==1)
				if (wavetype($s_test,1)==1)
					if (wave_index>0)
						insertPoints v_dim_wave,1,reference_waves
						insertPoints v_dim_wave,1,reference_waves_fit
					endif
					s_order+=s_test+";"
					s_nameofwave=testfolder+":reference"+num2str(wave_index)
					make/D/O w_inputwave
					killwaves w_inputwave
					wave w_inputwave=$s_test
					make/D/N=(v_rows_references)/O w_scaledwave
					for (index2=0;index2<v_rows_references;index2+=1)				//get the same scale
						v_position=energyscale[index2]
						v_position=((v_position - DimOffset(w_inputwave,0)) / DimDelta(w_inputwave,0))
						w_scaledwave[index2]=w_inputwave[v_position]
					endfor
					setscale/P x,v_sampleoffset,v_sampledelta,w_scaledwave
					duplicate w_scaledwave $s_nameofwave
					killwaves w_scaledwave
					reference_waves[wave_index]=$s_test
					reference_waves_fit[wave_index]=$s_nameofwave
					v_dim_wave=dimsize(reference_waves,0)
					wave_index+=1
				endif
			endif
		endfor
		lcm_variables[2]=wave_index
		v_num_of_waves=wave_index
		setdataFolder $testfolder
		if (wave_index<1)
			print "LCM: No reference waves."
			return 0	
		endif
	
	
	//transfer mask wave
	
	if (waveexists($s_mask_position)==1)
		wave w_mask_oldscale=$s_mask_position
		make/D/N=(v_rows_references)/O/FREE w_mask_newscale
		variable v_max_position=dimsize(w_mask_oldscale,0)-1
		for (index2=0;index2<v_rows_references;index2+=1)				//get the same scale
			v_position=energyscale[index2]
			v_position=((v_position - DimOffset(w_inputwave,0)) / DimDelta(w_inputwave,0))
			if ((v_position<0)||(v_position>v_max_position))
				w_mask_newscale[index2]=1
			else	
				w_mask_newscale[index2]=w_mask_oldscale[v_position]
			endif
		endfor
		setscale/P x,v_sampleoffset,v_sampledelta,w_mask_newscale		
	else
		v_use_mask=0
	endif
	
	//generate MT fit wave
		make/D/O/N=(v_num_of_waves,v_layer_sample) w_references_MT
		for (index=0;index<v_num_of_waves;index+=1)
			s_nameofwave=testfolder+":reference"+num2str(index)
			if (waveexists($s_nameofwave)==0)
				print "Error 4"
				print s_nameofwave
				return 0
			endif
			wave w_inputwave=$s_nameofwave
			for (index2=0;index2<v_layer_sample;index2+=1)
				w_references_MT[index][index2]=w_inputwave[index2]
			endfor
		endfor
		setscale/P y,v_sampleoffset,v_sampledelta,w_references_MT	
		
		
	//generate preview wave
		//print v_layer,v_columns,v_rows
		variable r_index, c_index, l_index
		if (v_chunks>1)
			make/N=(v_rows,(v_columns*v_layer),v_chunks)/FREE w_3dim
			variable ch_index
			for (r_index=0;r_index<v_rows;r_index++)
				for (c_index=0;c_index<v_columns;c_index++)
					for (l_index=0;l_index<v_layer;l_index++)
						for (ch_index=0;ch_index<v_chunks;ch_index++)
							w_3dim[r_index][(c_index*v_layer)+l_index][ch_index]=w_sample[r_index][c_index][l_index][ch_index]
						endfor
					endfor
				endfor	
			endfor
			make/D/O W_ISBeamAvg
			make/D/O W_ISBeamMax
			make/D/O W_ISBeamMin
			ImageStats/BEAM/Q w_3dim 
			setscale/P x,dimoffset(w_sample,3),dimdelta(w_sample,3),W_ISBeamAvg
			duplicate/O W_ISBeamAvg root:Packages:LCM:w_example
			killwaves W_ISBeamAvg
			killwaves W_ISBeamMax
			killwaves W_ISBeamMin			
		else	
			if (v_layer>0)
				if (v_columns>1&&v_rows>1)
					make/D/O W_ISBeamAvg
					make/D/O W_ISBeamMax
					make/D/O W_ISBeamMin
					ImageStats/BEAM/Q w_sample 
					setscale/P x,dimoffset(w_sample,2),dimdelta(w_sample,2),W_ISBeamAvg
					make/D/O root:Packages:LCM:w_example
					killwaves root:Packages:LCM:w_example
					duplicate W_ISBeamAvg root:Packages:LCM:w_example
					killwaves W_ISBeamAvg
					killwaves W_ISBeamMax
					killwaves W_ISBeamMin	
				else
					if (v_columns==1&&v_rows==1)
						make/D/O/N=(v_layer) w_build
						matrixOP/O w_build=beam(w_sample,0,0)
						setscale/P x,dimoffset(w_sample,2),dimdelta(w_sample,2),w_build
						duplicate/O w_build root:Packages:LCM:w_example								
					else
						if (v_rows>1)
							make/D/O/N=(v_layer) w_build
							matrixOP/O w_build=sumCols(w_sample)
							make/D/O/N=(v_layer) w_build2
							matrixOP/O w_build2=beam(w_build,0,0)
							setscale/P x,dimoffset(w_sample,2),dimdelta(w_sample,2),w_build2
							duplicate/O w_build2 root:Packages:LCM:w_example					
						else
							make/D/O/N=(v_layer) w_build
							matrixOP/O w_build=sumRows(w_sample)
							make/D/O/N=(v_layer) w_build2
							matrixOP/O w_build2=beam(w_build,0,0)
							setscale/P x,dimoffset(w_sample,2),dimdelta(w_sample,2),w_build2
							duplicate/O w_build2 root:Packages:LCM:w_example									
						endif
				
					endif

				endif		
			else
				if (v_columns>0)
					make/D/O/N=(v_columns) w_build
					matrixOP/O w_build=sumCols(w_sample)
					make/D/O/N=(v_columns) w_build2
					for (index=0;index<v_columns;index++)
						w_build2[index]=w_build[0][index]
					endfor
					setscale/P x,dimoffset(w_sample,1),dimdelta(w_sample,1),w_build2
					duplicate/O w_build2 root:Packages:LCM:w_example	
				else
					duplicate/O w_sample root:Packages:LCM:w_example			
				endif
			endif
		endif
	//correction for the case that start and end of fit is different	
		if (v_start_restart<v_sampleoffset)
			v_start_restart=v_sampleoffset
		endif
		if (v_end_restart<v_sampleoffset)
			v_end_restart=(v_sampleoffset+v_sampledelta*(v_rows_references-1))
		endif
		if (v_start_restart>(v_sampleoffset+v_sampledelta*(v_rows_references-1)))
			v_start_restart=v_sampleoffset
		endif
		if (v_end_restart>(v_sampleoffset+v_sampledelta*(v_rows_references-1)))
			v_end_restart=(v_sampleoffset+v_sampledelta*(v_rows_references-1))
		endif				
	//generate panel
		//lcm_variables //[0] = v_lcm_type: 1 = all-at-ones 2 = point-by-point , [1] = x_shift_done: 1 = no 2 = yes , [2] = quantity of waves, [3] = v_LCM_offset, [4] = v_LCM_delta, [5] = v_LCM_max, [6] = fitTol, [7] = start, [8] = end, [9] = dimensions of sample_wave
		make/D/O/N=(v_rows_references) root:Packages:LCM:w_mask=1
		make/D/O/N=(v_rows_references) root:Packages:LCM:w_mask_red=NaN
		wave w_mask=root:Packages:LCM:w_mask
		wave w_mask_red=root:Packages:LCM:w_mask_red
		if (v_use_mask>0)
			w_mask=w_mask_newscale
			for (index2=0;index2<dimsize(w_mask,0);index2+=1)
				if (numtype(w_mask_newscale[index2])!=0)
					w_mask_red[index2]=1
				endif
			endfor
		endif
		wave w_mask=root:Packages:LCM:w_mask
		wave w_mask_red=root:Packages:LCM:w_mask_red		
		setscale/P x,v_sampleoffset,v_sampledelta,w_mask 
		setscale/P x,v_sampleoffset,v_sampledelta,w_mask_red
		variable v_x_shift_abs_max=floor(((v_sampleoffset+v_sampledelta*(v_rows_references-1))-v_sampleoffset)*0.45)
		variable v_x_shift_start_max=floor(((v_sampleoffset+v_sampledelta*(v_rows_references-1))-v_sampleoffset)*0.2)
		lcm_variables[3]=v_sampleoffset
		lcm_variables[4]=v_sampledelta
		lcm_variables[5]=(v_sampleoffset+v_sampledelta*(v_rows_references-1))	
	//screen resolution	
		string s_igor_info=IgorInfo(0)
		variable v_info_elements=itemsinList(s_igor_info,";")
		string s_screen_info=replaceString("RECT=",stringfromList(v_info_elements-1,s_igor_info,";"),"")
		variable v_screen1_width=abs(str2num(stringfromList(3,s_screen_info,","))-str2num(stringfromList(1,s_screen_info,","))) //relative position
		variable v_screen1_length=abs(str2num(stringfromList(4,s_screen_info,","))-str2num(stringfromList(2,s_screen_info,",")))
		variable v_sc1_width_mp=round(v_screen1_width/2)
		variable v_sc1_length_mp=round(v_screen1_length/2)
		if (stringmatch(IgorInfo(2),"Windows")==1)	//if the OS is windows, the middle point of the Igor windows is used
			GetWindow kwFrameInner, wsize
			v_sc1_width_mp=round((V_right-V_left)/2) //relative position
			v_sc1_length_mp=round((V_bottom-V_top)/2)
		endif	
		if (stringmatch(IgorInfo(2),"Macintosh")==1) //if the OS is macOS, the middle point of the monitor that contains the mouse cursor is used
			GetMouse
			string s_Igorinfo = IgorInfo(0)
			string s_str
			variable v_n_screens=NumberByKey("NSCREENS", s_Igorinfo)
			//print v_n_screens
			variable v_i, v_screen_left, v_screen_top, v_screen_right, v_screen_bottom
			make/FREE/O/N=(v_n_screens,4) w_monitors
			for (v_i=1;v_i<=v_n_screens;v_i++)
				s_str = StringByKey("SCREEN"+num2str(v_i), s_Igorinfo)
				sscanf s_str, "DEPTH=%*d,RECT=%d,%d,%d,%d", v_screen_left, v_screen_top, v_screen_right, v_screen_bottom
				//print v_screen_left, v_screen_top, v_screen_right, v_screen_bottom
				w_monitors[v_i-1][0]=v_screen_left
				w_monitors[v_i-1][1]=v_screen_right
				w_monitors[v_i-1][2]=v_screen_top
				w_monitors[v_i-1][3]=v_screen_bottom
			endfor
			variable v_active_screen=1
			for (v_i=1;v_i<=v_n_screens;v_i++)
				if ((v_left>=w_monitors[v_i-1][0])&&(v_left<=w_monitors[v_i-1][1])&&(v_top>=w_monitors[v_i-1][2])&&(v_top<=w_monitors[v_i-1][3]))
					v_active_screen=v_i 
				endif
			endfor
			//print v_active_screen
			s_str = StringByKey("SCREEN"+num2str(v_active_screen), s_Igorinfo)
			sscanf s_str, "DEPTH=%*d,RECT=%d,%d,%d,%d", v_screen_left, v_screen_top, v_screen_right, v_screen_bottom //%*d ignors this iput, only %d is used
			//print v_screen_left, v_screen_right
			//print v_screen_top, v_screen_bottom
			v_sc1_width_mp=round((v_screen_right+v_screen_left)/2) //absolute position, negative values possible
			v_sc1_length_mp=round((v_screen_bottom+v_screen_top)/2)
		endif	
		//print v_sc1_width_mp,	v_sc1_length_mp
	//window for data input
		variable v_height=30*(wave_index)
		if (v_height < (75))
			v_height=(75)
		endif
		variable v_panel_left=v_sc1_width_mp-483 //165
		variable v_panel_rigth=v_sc1_width_mp+483 //(340+425+165+200+1)
		variable v_panel_top=v_sc1_length_mp-floor((v_height+110)/2) //140
		variable v_panel_bottom=v_sc1_length_mp+ceil((v_height+110)/2) //(v_height+110+140)
		make/O/N=4 root:Packages:LCM:w_panel_resolution
		wave w_panel_resolution=root:Packages:LCM:w_panel_resolution //left,right,top,bottom
		w_panel_resolution[0]=v_panel_left
		w_panel_resolution[1]=v_panel_rigth
		w_panel_resolution[2]=v_panel_top
		w_panel_resolution[3]=v_panel_bottom
		NewPanel/K=2/W=(v_panel_left,v_panel_top,v_panel_rigth,v_panel_bottom)/N=LCMpanel as "linear-combination modeling input" // /W=(165,140,(425+165),(30*(quantity_waves)+110+140))
		if (stringmatch(IgorInfo(2),"Macintosh")==1)
			MoveWindow v_panel_left,v_panel_top,-1,-1 //bug fix for negative positions on Macintosh
		endif
		ModifyPanel/W=LCMpanel fixedSize=1
		DrawRect 422,-2,(422+200),(v_height+110+140+2)
		DrawText 10	,25, "General settings"
		SetVariable V_FitTol	,pos={160	,10}	,value=_NUM:v_fittol_restart, limits={1e-10,0.01,0}, help={"Accuracy of the fit - between 0.01 and 1e-10."}
		DrawText 125	,25, "FitTol"
		SetVariable V_StartFit	,pos={270	,10}	,value=_NUM:v_start_restart, limits={v_sampleoffset,(v_sampleoffset+v_sampledelta*(v_rows_references-1)),0}, help={"Start value (energy) for fit."}
		DrawText 240	,25, "Start"
		SetVariable V_EndFit	,pos={365	,10}	,value=_NUM:v_end_restart, limits={v_sampleoffset,(v_sampleoffset+v_sampledelta*(v_rows_references-1)),0}, help={"End value (energy) for fit."}
		DrawText 340,25, "End"
		string s_LCM_sample_global=nameofwave(w_sample)
		string s_LCM_sample_local=ReplaceString("'",s_LCM_sample_global,"") //if wave starts with bad character
		string s_samplename
		variable v_start=v_sampleoffset
		variable v_end=(v_sampleoffset+v_sampledelta*(v_rows_references-1))
		if (strlen(s_LCM_sample_local)<38) //for long wavenames
			s_samplename="Reference waves for "+s_LCM_sample_local+""
		else
			s_samplename="Reference waves for "+s_LCM_sample_local[0,5]+" ... "+s_LCM_sample_local[(strlen(s_LCM_sample_local)-33),(strlen(s_LCM_sample_local)-1)]+""
		endif
		DrawText 10	,45, s_samplename
		//print v_allatones
		//print v_xshift_restart
		PopupMenu popup0 pos={430,5},size={50.00,14.00},value="All-at-ones Fit;Point-by-point-Fit",help={"The All-at-ones Fit is faster than the Point-by-point-Fit."},mode=v_allatones,disable=2 //,popmatch=s_allatones
		PopupMenu popup1 pos={430,25},size={50.00,14.00},value="No x-shift analysis;Multiple x-shift analysis;Single x-shift analysis",help={"Better results for x-shift analysis if Initial Guess of waves >> 1e-6 (e. g. 0.5), ± = zero holds the initial value, Single x-shift analysis uses the initial x-shift of the first wave"},mode=v_xshift_restart //popmatch=s_xshift_restart,
		string s_start
		string s_min
		string s_max, s_x_start, s_x_var, s_x_hold
		string s_wavepos=""
		variable v_intemsinlist
		for (index=0;index<wave_index;index+=1)
			s_start="V_Start_"+num2str(index)
			s_min="V_Min_"+num2str(index)
			s_max="V_Max_"+num2str(index)
			s_x_start="V_xStart_"+num2str(index)
			s_x_var="V_xVar_"+num2str(index)
			s_wavepos=nameofwave(reference_waves[index])
			s_wavepos=ReplaceString("'",s_wavepos,"") //if wave starts with bad character
			DrawText 10	,((50+(30*index))+15), s_wavepos[0,12]
			SetVariable $s_start	,pos={175	,(50+(30*index))}	,value=_NUM:w_initial_guess[index], limits={-inf,inf,0} , help={"Initial Guess for the fit, must be non-zero, e.g. 1e-6, of the wave "+getwavesDataFolder(reference_waves[index],2)+"."}
			DrawText 105	,((50+(30*index))+15), "Initial Guess"
			SetVariable $s_min	,pos={270	,(50+(30*index))}	,value=_NUM:w_min[index], limits={-inf,inf,0} , help={"Lower limit of the fit, usually -inf or zero"}
			DrawText 242	,((50+(30*index))+15), "Min"
			SetVariable $s_max	,pos={365	,(50+(30*index))}	,value=_NUM:w_max[index], limits={-inf,inf,0} , help={"Upper limit of the fit, usually +inf"}
			DrawText 340,((50+(30*index))+15), "Max"
			SetVariable $s_x_start	,pos={500	,(50+(30*index))}	,value=_NUM:w_initial_xshift[index], limits={-v_x_shift_start_max,v_x_shift_start_max,0} , help={"Initial Guess, can be zero"}
			DrawText 430	,((50+(30*index))+15), "Initial x-shift"
			SetVariable $s_x_var	,pos={565	,(50+(30*index))}	,value=_NUM:w_var_xshift[index], limits={-v_x_shift_abs_max,v_x_shift_abs_max,0} , help={"Range for energyshift, e.g. "+num2str(10*v_sampledelta)+", Zero is interpreted as hold value"}
			DrawText 553	,((50+(30*index))+15), "±"		
		endfor
		
		DrawText 10	,((50+(30*index))+15), "absolute offset"
		SetVariable V_start_abs	,pos={175	,(50+(30*index))}	,value=_NUM:w_initial_guess[index], limits={-inf,inf,0}, help={"Initial Guess of the offset for the fit, e.g. 1e-6"}
		DrawText 105	,((50+(30*index))+15), "Initial Guess"
		SetVariable V_min_abs	,pos={270	,(50+(30*index))}	,value=_NUM:w_min[index], limits={-inf,inf,0}, help={"Lower limit of the offset, usually -inf"}
		DrawText 242	,((50+(30*index))+15), "Min"
		SetVariable V_max_abs	,pos={365	,(50+(30*index))}	,value=_NUM:w_max[index], limits={-inf,inf,0}, help={"Upper limit of the offset, usually +inf"}
		DrawText 340,((50+(30*index))+15), "Max"
		Button B_Doit			,pos={253	,(v_height)+80}	,size={75,20}		,title="Do It" ,fColor=(257*30,257*144,257*255)		,proc=f_Button_doit_LCM, help={"Performing the fit can take several minutes. The fitted wave is "+getwavesDataFolder(w_sample,2)+"."}
		Button B_Cancel			,pos={338	,(v_height)+80}	,size={75,20}		,title="Cancel"		,proc=f_Button_cancel_LCM
		Button B_SpectraPrev,pos={10	,(v_height)+80}	,size={55,20}		,title="Waves", fColor=(52428,52428,52428)	,proc= f_Button_wave_preview_LCM, help={"Preview of the used reference waves."}
		Button B_SaveConfig,pos={75	,(v_height)+80}	,size={100,20}		,title="Save Config.",proc = f_Button_saveconfig_LCM, help={"Save the visible configuration in root:."}
		DrawText 430	,(v_height)+80+18), "Set initial Guess"
		Button B_initial_a,pos={530	,(v_height)+80}	,size={40,20} ,title="1e-6" ,proc=f_Button_initial_1em6_LCM
		Button B_initial_b,pos={575	,(v_height)+80}	,size={40,20} ,title="0.5" ,proc=f_Button_initial_0_5_LCM
		DrawText 430	,(v_height)+55+18), "Min"
		Button B_min_a,pos={452	,(v_height)+55}	,size={40,20} ,title="-inf" ,proc=f_Button_min_minf_LCM
		Button B_min_b,pos={497	,(v_height)+55}	,size={40,20} ,title="0" ,proc=f_Button_min_0_LCM	
		DrawText 548	,(v_height)+55+18), "Max"
		Button B_max_b,pos={575	,(v_height)+55}	,size={40,20} ,title="+inf" ,proc=f_Button_max_inf_LCM	
	//window part for mask
		//SetDrawLayer UserBack
		Button B_mask,pos={422+200+10.00,155.00},size={100.00,20.00},title="Mask region",proc=f_Button_Mask_NaN_LCM, help={"Spectral region does not include in the fit (red region)"}
		Button B_mask,fColor=(65535,0,0)
		Button B_unmask,pos={422+200+120.00,155.00},size={100.00,20.00},title="Unmask region",proc=f_Button_Mask_1_LCM, help={"Spectral region does include in the fit (blue region)"}
		Button B_unmask,fColor=(0,43690,65535)
		SetVariable V_B_start,pos={422+200+30.00,129.00},size={100.00,14.00},title="Region Start",value=_NUM:v_start, help={"For individual points, the start and end values must be set equal."}
		SetVariable V_B_start,limits={v_start,v_end,0}
		SetVariable V_B_end,pos={422+200+165.00,129.00},size={100.00,14.00},title="Region End",value=_NUM:v_end
		SetVariable V_B_end,limits={v_start,v_end,0}
		Button B_unmask_all,pos={422+200+230.00,155.00},size={100.00,20.00},title="Unmask all", proc=f_Button_Mask_reset_LCM, help={"The mask is reset so that all spectral regions are included in the fit."}
		Button B_unmask_all,fColor=(1,65535,33232)
		Display/L/W=(422+200+10,9,422+200+10+323,119)/HOST=# w_mask
		AppendToGraph/L w_mask_red
		AppendToGraph/R root:Packages:LCM:w_example			
		ModifyGraph hbFill=2
		ModifyGraph lSize=1.5
		SetAxis left 0,1
		ModifyGraph mode(w_mask)=7,lsize(w_mask)=0,plusRGB(w_mask)=(1,16019,65535,32768),negRGB(w_mask)=(1,16019,65535,32768)
		ModifyGraph mode(w_mask_red)=7,lsize(w_mask_red)=0,plusRGB(w_mask_red)=(65535,0,0,32768),negRGB(w_mask_red)=(65535,0,0,3)
		ModifyGraph hbFill(w_mask)=3	
		ModifyGraph rgb(w_example)=(0,0,0)
		ModifyGraph usePlusRGB=1
		ModifyGraph useNegRGB=1
		ModifyGraph mirror(bottom)=2
		ModifyGraph notation(left)=1
		ModifyGraph nticks(left)=0
		ModifyGraph noLabel(left)=1
		ModifyGraph nticks(right)=0
		ModifyGraph noLabel(right)=1
		ModifyGraph UIControl=(2^0+2^1+2^3+2^4+2^5+2^6+2^7+2^8+2^9+2^10+2^11)
		Label bottom ""
		Label left "Average"				
		ModifyGraph margin(right)=15,margin(top)=15,margin(left)=20
		RenameWindow #,G0
		SetActiveSubwindow ##		
end

//Button to cancel the function
Function f_Button_cancel_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba	
	switch( ba.eventCode )
		case 2: // mouse up
				if (WinType("maskinput"))
					DoWindow/F maskinput
					return 0
				else
					if (WinType("spectraoverview"))
						killwindow spectraoverview
					endif	
					print "LCM: Abort by user."					
					killwindow/Z LCMpanel
					setdataFolder root:
					killdataFolder/Z root:Packages:LCM:
				endif
				break
		case -1: // control being killed
			break
	endswitch
	return 0
end

//Buttons to change multiple parameters of the input
Function f_Button_initial_1em6_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	switch( ba.eventCode )
		case 2: // mouse up				
			f_change_LCM_start_values(1e-6)
			break
		case -1: // control being killed
			break
	endswitch
	return 0
End

Function f_Button_initial_0_5_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	switch( ba.eventCode )
		case 2: // mouse up				
			f_change_LCM_start_values(0.5)
			break
		case -1: // control being killed
			break
	endswitch
	return 0
End

Function f_Button_min_minf_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	switch( ba.eventCode )
		case 2: // mouse up				
			f_change_LCM_min_values_all(-inf)
			break
		case -1: // control being killed
			break
	endswitch
	return 0
End

Function f_Button_min_0_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	switch( ba.eventCode )
		case 2: // mouse up				
			f_change_LCM_min_values(0)
			break
		case -1: // control being killed
			break
	endswitch
	return 0
End

Function f_Button_max_inf_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	switch( ba.eventCode )
		case 2: // mouse up				
			f_change_LCM_max_values(inf)
			break
		case -1: // control being killed
			break
	endswitch
	return 0
End

Function f_Button_max_0_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	switch( ba.eventCode )
		case 2: // mouse up				
			f_change_LCM_max_values(0)
			break
		case -1: // control being killed
			break
	endswitch
	return 0
End

//Buttons for the wave prewview, blue = wave that will be fitted, black = reference waves
function f_Button_wave_preview_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	wave/WAVE w_reference_waves_fit = root:Packages:LCM:reference_waves_fit
	wave/WAVE w_reference_waves = root:Packages:LCM:reference_waves
	wave w_example = root:Packages:LCM:w_example
	variable v_waves=dimsize(root:Packages:LCM:reference_waves_fit,0)+1
	variable index, v_start, v_end, v_scale_start, v_scale_end
	string axis_name
	switch( ba.eventCode )
		case 2: // mouse up
			if (WinType("spectraoverview"))
				DoWindow/F spectraoverview
				return 0
			else
				variable inputvar
				wave w_panel_resolution=root:Packages:LCM:w_panel_resolution //left,right,top,bottom
				variable v_mainpanel_rigth=w_panel_resolution[1]
				variable v_mainpanel_top=w_panel_resolution[2]
				variable v_distance=2
				if (stringmatch(IgorInfo(2),"Windows")==1)
					v_distance=8
					//print v_distance
				endif
				NewPanel /K=2 /W=((v_mainpanel_rigth+v_distance),(v_mainpanel_top),(v_mainpanel_rigth+v_distance+340),(v_mainpanel_top+560)) /N=spectraoverview as "linear-combination reference waves" //(340+425+165+200+2),(140),(1472),(700)
				ModifyPanel /W=spectraoverview ,noEdit=1 ,fixedSize=1  
				SetDrawLayer UserBack
				Button B_done,pos={230.00,530.00},size={100.00,20.00},title="Close",proc=f_Button_spectra_done_LCM
				Button B_done,fColor=(52428,52428,52428)
				Display/L/W=(2,9,361,250+270)/HOST=#  w_example
				ModifyGraph axisEnab(left)={1-(1/v_waves),1}
				Label left "Average"
				ModifyGraph nticks(left)=0
				ModifyGraph noLabel(left)=1
				v_scale_start=dimOffset(w_example,0)
				v_scale_end=dimOffset(w_example,0)+dimdelta(w_example,0)*(dimsize(w_example,0)-1)
				for (index=0;index<(v_waves-1);index++)
					axis_name="axis"+num2str(index)
					//print axis_name
					appendtograph/L=$axis_name w_reference_waves_fit[index]
					v_start=(1-(index+2)*(1/v_waves))
					//print v_start
					v_end=(1-(index+1)*(1/v_waves))
					//print v_end
					ModifyGraph nticks($axis_name)=0
					ModifyGraph noLabel($axis_name)=1
					ModifyGraph axisEnab($axis_name)={v_start,v_end}
					ModifyGraph freePos($axis_name)=0
					ModifyGraph lblPosMode($axis_name)=2;DelayUpdate
					Label $axis_name nameofwave(w_reference_waves[index])
				endfor	
				SetAxis bottom v_scale_start,v_scale_end
				ModifyGraph lSize=1.5
				ModifyGraph hbFill=2
				ModifyGraph mirror(bottom)=1
				ModifyGraph notation(left)=1
				ModifyGraph margin(left)=20
				ModifyGraph rgb=(0,0,0)
				ModifyGraph rgb(w_example)=(1,16019,65535)
				Label bottom ""
				ModifyGraph mirror=0
				ModifyGraph tick(left)=0
				ModifyGraph tick(bottom)=0
				NewFreeAxis/O/R mirrorR
				ModifyFreeAxis mirrorR, master=left
				ModifyGraph freePos(mirrorR)={0,kwFraction}
				ModifyGraph tick(mirrorR)=2,noLabel(mirrorR)=2
				ModifyGraph nticks(mirrorR)=0
				NewFreeAxis/O/T mirrorT
				ModifyFreeAxis mirrorT, master=bottom
				ModifyGraph freePos(mirrorT)={0,kwFraction}
				ModifyGraph tick(mirrorT)=2,noLabel(mirrorT)=2
				ModifyGraph margin(right)=15,margin(top)=15			
				RenameWindow #,G1
				SetActiveSubwindow ##
			endif
			break
		case -1: // control being killed
			break
	endswitch
end

Function f_Button_spectra_done_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	switch( ba.eventCode )
		case 2: // mouse up				
			killwindow spectraoverview
			break
		case -1: // control being killed
			break
	endswitch
	return 0
End

//Button to unmask the complete region / mask_wave
Function f_Button_Mask_reset_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	wave w_mask=root:Packages:LCM:w_mask
	wave w_mask_red=root:Packages:LCM:w_mask_red
	switch( ba.eventCode )
		case 2: // mouse up
			w_mask=1
			w_mask_red=NaN
			break
		case -1: // control being killed
			break
	endswitch
	return 0
End

//Button to mask an region
Function f_Button_Mask_NaN_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	wave w_mask=root:Packages:LCM:w_mask
	wave w_mask_red=root:Packages:LCM:w_mask_red
	variable v_rows=dimsize(w_mask,0)-1
	variable v_start
	variable v_end
	variable b_index
	switch( ba.eventCode )
		case 2: // mouse up
			ControlInfo/W=LCMpanel V_B_start
			v_start=ScaleToIndex(w_mask,V_value,0)
			ControlInfo/W=LCMpanel V_B_end
			v_end=ScaleToIndex(w_mask,V_value,0)
			if (v_start>v_end)
				b_index=v_start
				v_start=v_end
				v_end=b_index
			endif
			if (v_start<0)
				v_start=0
			endif
			if (v_end<0)
				v_end=0
			endif
			if (v_start>v_rows)
				v_start=v_rows
			endif	
			if (v_end>v_rows)
				v_end=v_rows
			endif	
			if (v_start==v_end)
				w_mask[v_start]=NaN
				w_mask_red[v_start]=1
			else
				for (b_index=v_start;b_index<(v_end+1);b_index+=1)
					w_mask[b_index]=NaN
					w_mask_red[b_index]=1
				endfor
			endif
			break
		case -1: // control being killed
			break
	endswitch
	return 0
End

//Button to unmask an region
Function f_Button_Mask_1_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	wave w_mask=root:Packages:LCM:w_mask
	wave w_mask_red=root:Packages:LCM:w_mask_red
	variable v_rows=dimsize(w_mask,0)-1
	variable v_start
	variable v_end
	variable b_index
	switch( ba.eventCode )
		case 2: // mouse up
			ControlInfo/W=LCMpanel V_B_start
			v_start=ScaleToIndex(w_mask,V_value,0)
			ControlInfo/W=LCMpanel V_B_end
			v_end=ScaleToIndex(w_mask,V_value,0)
			if (v_start>v_end)
				b_index=v_start
				v_start=v_end
				v_end=b_index
			endif
			if (v_start<0)
				v_start=0
			endif
			if (v_end<0)
				v_end=0
			endif
			if (v_start>v_rows)
				v_start=v_rows
			endif	
			if (v_end>v_rows)
				v_end=v_rows
			endif						
			if (v_start==v_end)
				w_mask[v_start]=1
				w_mask_red[v_start]=NaN
			else
				for (b_index=v_start;b_index<(v_end+1);b_index+=1)
					w_mask[b_index]=1
					w_mask_red[b_index]=NaN
				endfor
			endif
			break
		case -1: // control being killed
			break
	endswitch
	return 0
End

function f_LCM_save_config()
	variable v_num_of_waves=dimsize(root:Packages:LCM:reference_waves,0)
 	variable v_information_rows
 	if ((v_num_of_waves+3)>9)
		v_information_rows=v_num_of_waves+3
	else
		v_information_rows=9
	endif	
 	
 	variable V_value	
	ControlInfo/W=LCMpanel popup1 //1 = no x-shift, 2 = multiple x-shift, 3 = single x-shift
	variable x_shift_done=V_value
	switch (x_shift_done)
		case 1:
			make/FREE/T/O/N=((v_information_rows),5) w_config
		break
		case 2:
			make/FREE/T/O/N=((v_information_rows),7) w_config
	 		w_config[0][5]="Initial x-shift"
	 		w_config[0][6]="±"			
		break
		case 3:
			make/FREE/T/O/N=((v_information_rows),7) w_config
	 		w_config[0][5]="Initial x-shift"
	 		w_config[0][6]="±"		
		break
	endswitch	
 	wave/WAVE reference_waves=root:Packages:LCM:reference_waves
 	variable v_fitTol, v_start, v_end, v_start_energy, v_end_energy, index
 	string s_fit_type,s_ref_waves_names

	ControlInfo/W=LCMpanel popup0 //1 = all-at-ones, 2 = point-by-point
	switch (V_value)
		case 1:
			s_fit_type="All-at-ones Fit"		
		break
		case 2:
			s_fit_type="Point-by-point Fit"		//don't use, only for testing of new funcfit functions (Point-by-point forgives more mistakes)
		break
	endswitch
	ControlInfo/W=LCMpanel V_FitTol
	v_fitTol=V_value
	ControlInfo/W=LCMpanel V_StartFit
	v_start=V_value
	ControlInfo/W=LCMpanel V_EndFit
	v_end=V_value
 	w_config[0][0]="0"
	w_config[0][0]=Secs2Date(DateTime,2)
	w_config[1][0]=Secs2Time(DateTime,2)
	w_config[2][0]="Saved configuration"
	w_config[6][0]="FitTol: "+num2str(v_fitTol)
	w_config[5][0]="Scale: "+num2str(v_start)+" to "+num2str(v_end)
	w_config[8][0]=s_fit_type
	w_config[0][1]="References"
	w_config[0][2]="Initial Guess"
	w_config[0][3]="Min"
	w_config[0][4]="Max" 	
	for (index=0;index<(v_num_of_waves);index++)
		w_config[index+1][1]=GetWavesDataFolder(reference_waves[index],2)
		s_ref_waves_names+=nameofwave(reference_waves[index])+";"
	endfor
	string s_transfer
	for (index=0;index<(v_num_of_waves);index++)
		//print index
		s_transfer="V_Start_"+num2str(index)
		//print s_transfer
		ControlInfo/W=LCMpanel $s_transfer
		w_config[index+1][2]=num2str(V_Value)
		s_transfer="V_Min_"+num2str(index)
		//print s_transfer
		ControlInfo/W=LCMpanel $s_transfer
		w_config[index+1][3]=num2str(V_Value)
		s_transfer="V_Max_"+num2str(index)
		//print s_transfer
		ControlInfo/W=LCMpanel $s_transfer
		w_config[index+1][4]=num2str(V_Value)
		if (x_shift_done>1)
			s_transfer="V_xStart_"+num2str(index)
			//print s_transfer
			ControlInfo/W=LCMpanel $s_transfer
			w_config[index+1][5]=num2str(V_Value)
			s_transfer="V_xVar_"+num2str(index)
			//print s_transfer
			ControlInfo/W=LCMpanel $s_transfer
			w_config[index+1][6]=num2str(V_Value)
		endif
	endfor	
		//print index
		ControlInfo/W=LCMpanel V_Start_abs
		w_config[v_num_of_waves+1][2]=num2str(V_Value)
		ControlInfo/W=LCMpanel V_Min_abs
		w_config[v_num_of_waves+1][3]=num2str(V_Value)
		ControlInfo/W=LCMpanel V_Max_abs	
		w_config[v_num_of_waves+1][4]=num2str(V_Value)	
		w_config[v_num_of_waves+1][1]="absolute Offset"
		s_ref_waves_names+="abs;"
	switch (x_shift_done)
		case 1:
		break
		case 2:
	 		w_config[8][0]+=" + Multiple x-shift analysis"			
		break
		case 3:
	 		w_config[8][0]+=" + Single x-shift analysis"			
		break
	endswitch	
	string s_export
	for (index=1;index<inf;index++)
		s_export="root:'"+Secs2Date(DateTime,-2,"_")+"_"+num2str(index)+"_LCM_configuration'"
		if (waveexists($s_export)==0)
			index=inf
		endif
	endfor
	duplicate/O w_config $s_export
end

function f_Button_saveconfig_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	switch( ba.eventCode )
		case 2: // mouse up
			f_LCM_save_config()
		break
		case -1: // control being killed
			break
	endswitch
	return 0
End

Function f_Button_doit_LCM(ba) : ButtonControl
	STRUCT WMButtonAction &ba
	wave/WAVE w_ref_waves=root:Packages:LCM:reference_waves
	variable index
	string s_control
	switch( ba.eventCode )
		case 2: // mouse up
			if (WinType("spectraoverview"))
				killwindow spectraoverview
			endif
			variable v_waves_num=dimsize(w_ref_waves,0)
			make/D/O/N=(v_waves_num+1) root:Packages:LCM:min_wave
			wave min_wave=root:Packages:LCM:min_wave
			make/D/O/N=(v_waves_num+1) root:Packages:LCM:max_wave
			wave max_wave=root:Packages:LCM:max_wave
			make/D/O/N=(v_waves_num) root:Packages:LCM:var_x_wave
			wave var_x_wave=root:Packages:LCM:var_x_wave
			make/D/O/N=(v_waves_num+1) root:Packages:LCM:start_val
			wave start_val=root:Packages:LCM:start_val
			make/D/O/N=(v_waves_num) root:Packages:LCM:start_x_val
			wave start_x_val=root:Packages:LCM:start_x_val
			//s_control=ControlNameList("LCMpanel", ";", "V*")
			wave lcm_variables=root:Packages:LCM:lcm_variables //[0] = v_lcm_type: 1 = all-at-ones 2 = point-by-point , [1] = x_shift_done: 1 = no 2 = yes , [2] = quantity of waves, [3] = v_LCM_offset, [4] = v_LCM_delta, [5] = v_LCM_max, [6] = fitTol, [7] = start, [8] = end , [9] = dimensions of sample_wave
			ControlInfo/W=LCMpanel popup0 //1 = all-at-ones, 2 = point-by-point
			lcm_variables[0]=V_value
			ControlInfo/W=LCMpanel popup1 //1 = no x-shift, 2 = multiple x-shift, 3 = single x-shift
			lcm_variables[1]=V_value
			ControlInfo/W=LCMpanel V_FitTol
			lcm_variables[6]=V_value
			ControlInfo/W=LCMpanel V_StartFit
			lcm_variables[7]=V_value
			ControlInfo/W=LCMpanel V_EndFit
			lcm_variables[8]=V_value
			string s_transfer		
			for (index=0;index<((dimsize(w_ref_waves,0)));index++)
				//print index
				s_transfer="V_Start_"+num2str(index)
				//print s_transfer
				ControlInfo/W=LCMpanel $s_transfer
				start_val[index]=V_Value
				s_transfer="V_Min_"+num2str(index)
				//print s_transfer
				ControlInfo/W=LCMpanel $s_transfer
				min_wave[index]=V_Value
				s_transfer="V_Max_"+num2str(index)
				//print s_transfer
				ControlInfo/W=LCMpanel $s_transfer
				max_wave[index]=V_Value
				s_transfer="V_xStart_"+num2str(index)
				//print s_transfer
				ControlInfo/W=LCMpanel $s_transfer
				start_x_val[index]=V_Value
				s_transfer="V_xVar_"+num2str(index)
				//print s_transfer
				ControlInfo/W=LCMpanel $s_transfer
				var_x_wave[index]=V_Value
			endfor
			//V_Start_0;V_Min_0;V_Max_0;V_xStart_0;V_xVar_0;V_start_abs;V_min_abs;V_max_abs
			ControlInfo/W=LCMpanel V_start_abs
			start_val[index]=V_Value
			ControlInfo/W=LCMpanel V_min_abs
			min_wave[index]=V_Value
			ControlInfo/W=LCMpanel V_max_abs
			max_wave[index]=V_Value
			killwindow LCMpanel
			if (waveexists(root:Packages:LCM:reference_waves)==0)
				print "Error, reference wave(s) missing."
				return 0
			endif
			if (waveexists(root:Packages:LCM:lcm_variables)==0)
				print "Error, variables wave missing."
				return 0
			endif
			if (waveexists(root:Packages:LCM:w_references_MT)==0)
				print "Error, MD reference wave(s) missing."
				return 0
			endif			
			f_do_fit_LCM()
			break
		case -1: // control being killed
			break
	endswitch
	return 0
End

function f_do_fit_LCM()
	setdatafolder root:Packages:LCM:
	//input of all temp. waves
	wave/WAVE w_ref_waves=root:Packages:LCM:reference_waves
	wave start_val=root:Packages:LCM:start_val
	wave min_wave=root:Packages:LCM:min_wave
	wave max_wave=root:Packages:LCM:max_wave
	wave start_x_val=root:Packages:LCM:start_x_val
	wave var_x_wave=root:Packages:LCM:var_x_wave
	wave/WAVE w_data_wave=root:Packages:LCM:data_wave
	wave w_sample_input=w_data_wave[0]
	wave/WAVE reference_waves=root:Packages:LCM:reference_waves
	wave/WAVE reference_waves_fit=root:Packages:LCM:reference_waves_fit
	wave w_references_MT=root:Packages:LCM:w_references_MT
	wave w_mask=root:Packages:LCM:w_mask
	wave lcm_variables=root:Packages:LCM:lcm_variables //[0] = v_lcm_type: 1 = all-at-ones 2 = point-by-point , [1] = x_shift_done: 1 = no 2 = yes (multiple) 3 = yes (single) , [2] = quantity of waves, [3] = v_LCM_offset, [4] = v_LCM_delta, [5] = v_LCM_max, [6] = fitTol, [7] = start, [8] = end , [9] = dimensions of sample_wave
	//read of input_variables
	variable v_lcm_type_inp=lcm_variables[0] //1 = all-at-ones 2 = point-by-point
	//print v_lcm_type_inp
	variable x_shift_done=lcm_variables[1] //1 = no 2 = yes (multiple) 3 = yes (single) 
	//print x_shift_done
	variable v_num_of_waves=lcm_variables[2]
	variable v_LCM_offset=lcm_variables[3]
	variable v_LCM_delta=lcm_variables[4]
	variable v_LCM_max=lcm_variables[5]
	variable v_fitTol=lcm_variables[6]
	variable v_start=lcm_variables[7]
	variable v_end=lcm_variables[8]
	variable v_inputwave_dim=lcm_variables[9]
	//generate needed variables
	variable v_threads						//How many threads?
	variable v_rows_max						//How many rows?
	variable v_columns_max
	variable v_layers_max
	variable v_chunk_max
	variable v_layers_old=NaN
	variable v_row_per_thread		//How many rows for each thread
	variable r_index,c_index,l_index,t_index,index,ch_index	
	variable v_rows, v_columns, v_layers
	//create 3-dim. wave for 1-dim. and 2-dim inputs, thread group parameters
	switch(v_inputwave_dim)
		case 1: //1D (spectrum = rows)
			make/D/O/N=(1,1,dimsize(w_sample_input,0)) w_sample=0
			for (index=0;index<(dimsize(w_sample_input,0));index++)
				w_sample[0][0][index]=w_sample_input[index]
			endfor
			setscale/P z,dimoffset(w_sample_input,0),dimdelta(w_sample_input,0),w_sample
		 	v_threads=1						//How many threads? -> 1-dim. does not need MT
			v_rows_max=dimsize(w_sample,0)						//How many rows?
		 	v_columns_max=dimsize(w_sample,1)
		 	v_layers_max=dimsize(w_sample,2)
		 	v_row_per_thread=v_rows_max			//How many rows for each thread -> 1-dim.: one thread does everything						
		break
		
		case 2:	//2D (spectra = rows, positions = columns)
			v_rows=dimsize(w_sample_input,0)
			v_columns=dimsize(w_sample_input,1)		
			make/D/N=(v_rows,1,v_columns) w_sample=0
			for (l_index=0;l_index<v_columns;l_index++)
				for (r_index=0;r_index<v_rows;r_index++)
					w_sample[r_index][0][l_index]=w_sample_input[r_index][l_index]
				endfor
			endfor
			setscale/P z,dimoffset(w_sample_input,1),dimdelta(w_sample_input,1),w_sample
		 	v_threads=ThreadProcessorCount						//How many threads?
			v_rows_max=dimsize(w_sample,0)						//How many rows?
		 	v_columns_max=dimsize(w_sample,1)
		 	v_layers_max=dimsize(w_sample,2)
		 	v_row_per_thread=ceil(v_rows_max/v_threads)			//How many rows for each thread			
		break
		
		case 3:	//3D
			duplicate w_sample_input w_sample
		 	v_threads=ThreadProcessorCount						//How many threads?
			v_rows_max=dimsize(w_sample,0)						//How many rows?
		 	v_columns_max=dimsize(w_sample,1)
		 	v_layers_max=dimsize(w_sample,2)
		 	v_row_per_thread=ceil(v_rows_max/v_threads)			//How many rows for each thread			
		break
		
		case 4: //4D
			v_rows_max=dimsize(w_sample_input,0)
			v_columns_max=dimsize(w_sample_input,1)
			v_layers_max=dimsize(w_sample_input,2)
			v_chunk_max=dimsize(w_sample_input,3)
			v_layers_old=dimsize(w_sample_input,2)
			make/D/N=(v_rows_max,v_columns_max*v_layers_old,v_chunk_max)/O w_sample
			for (r_index=0;r_index<v_rows_max;r_index++)
				for (c_index=0;c_index<v_columns_max;c_index++)
					for (l_index=0;l_index<v_layers_old;l_index++)
						for (ch_index=0;ch_index<v_chunk_max;ch_index++)
							w_sample[r_index][(c_index*v_layers_old)+l_index][ch_index]=w_sample_input[r_index][c_index][l_index][ch_index]
						endfor
					endfor
				endfor	
			endfor
			v_layers_max=v_chunk_max
			//duplicate/O w_sample root:w_sample_test
			v_rows_max=dimsize(w_sample,0)
			v_columns_max=dimsize(w_sample,1)
			v_layers_max=dimsize(w_sample,2)
			v_threads=ThreadProcessorCount
			v_row_per_thread=ceil(v_rows_max/v_threads)
		break
	endswitch
	//duplicate/O w_sample root:w_sample_test
	v_rows=dimsize(w_sample,0)
	v_columns=dimsize(w_sample,1)
	v_layers=dimsize(w_sample,2)
	//print v_rows,v_columns,v_layers
	//create export-waves & correction of variable v_lcm_type
	make/D/O/N=(v_rows,v_columns,v_layers) w_fits=NaN
	string s_fit_type
	variable v_lcm_type
	switch(v_lcm_type_inp)
		case 1:
			v_lcm_type=1
			s_fit_type="All-at-ones Fit"
		break
		case 2:
			v_lcm_type=0
			s_fit_type="Point-by-point Fit"
		break
	endswitch
	switch (x_shift_done)
		case 1: //no energy shift
			make/D/O/N=(v_rows,v_columns,(v_num_of_waves+1)) w_export=NaN
			make/D/O/N=(v_rows,v_columns,(v_num_of_waves+1)) w_stddev=NaN	
		break
		case 2: //energy shift for each wave individual
			v_lcm_type+=2
			make/D/O/N=(v_rows,v_columns,(2*v_num_of_waves+1)) w_export=NaN
			make/D/O/N=(v_rows,v_columns,(2*v_num_of_waves+1)) w_stddev=NaN	
		break
		case 3: //one energy shift for each wave together
			v_lcm_type+=4
			make/D/O/N=(v_rows,v_columns,(v_num_of_waves+2)) w_export=NaN
			make/D/O/N=(v_rows,v_columns,(v_num_of_waves+2)) w_stddev=NaN		
		break
	endswitch	
	//print "v_lcm_type "+num2str(v_lcm_type)
	//hold string
	string s_hold_LCM=""
	variable v_value_test=0
	variable v_xshift_test=0
	switch (x_shift_done)
		case 1: //no x-shift analysis
			for (index=0;index<(v_num_of_waves+1);index++)
				if (min_wave[index]==max_wave[index])
					s_hold_LCM+="1"
					v_value_test+=1
				else
					s_hold_LCM+="0"
				endif
			endfor
		break
		case 2: //x-shift analysis
			for (index=0;index<(v_num_of_waves+1);index++)	
				if (min_wave[index]==max_wave[index])
					s_hold_LCM+="1"
					v_value_test+=1
				else
					s_hold_LCM+="0"
				endif
			endfor	
			for (index=0;index<(v_num_of_waves);index++)	
				if (var_x_wave[index]==0)
					s_hold_LCM+="1"
					v_value_test+=1
				else
					s_hold_LCM+="0"
					v_xshift_test+=1
				endif
			endfor
			if (v_xshift_test==0)
				x_shift_done=0
			endif						
		break
		case 3:
			for (index=0;index<(v_num_of_waves+1);index++)	
				if (min_wave[index]==max_wave[index])
					s_hold_LCM+="1"
					v_value_test+=1
				else
					s_hold_LCM+="0"
				endif
			endfor
			s_hold_LCM+="0"	 //single x shift = never holded				
		break
	endswitch
	//print s_hold_LCM
	//create contr string (K values) + threadsafe transformation
	variable var_min, var_max, var_test	
	string s_constraint="\""
		for (index=0;index<(v_num_of_waves+1);index+=1)
			var_min=min_wave[index]
			var_max=max_wave[index]
			if (var_min>-inf&&var_min!=var_max)
				s_constraint+="K"+num2str(index)+" > "+num2str(var_min)+","
				var_test+=1
			endif
			if (var_max<inf&&var_min!=var_max)
				s_constraint+="K"+num2str(index)+" < "+num2str(var_max)+","
				var_test+=1
			endif
		endfor	
	 
	 if (x_shift_done==2) //if x-shift is done
		for (index=0;index<(v_num_of_waves);index+=1)
			var_min=-abs(var_x_wave[index])
			var_max=abs(var_x_wave[index])
			if (var_max!=0)
				if (var_min>-inf&&var_min!=var_max)
					s_constraint+="K"+num2str(v_num_of_waves+1+index)+" > "+num2str(var_min)+","
					var_test+=1
				endif
				if (var_max<inf&&var_min!=var_max)
					s_constraint+="K"+num2str(v_num_of_waves+1+index)+" < "+num2str(var_max)+","
					var_test+=1
				endif	
			endif	
		endfor
	endif	
	if (x_shift_done==3) //if x-shift is done
		var_min=-abs(var_x_wave[0])
		var_max=abs(var_x_wave[0])		 
		var_max=abs(var_x_wave[0])
		if (var_max!=0)
			if (var_min>-inf&&var_min!=var_max)
				s_constraint+="K"+num2str(v_num_of_waves+1)+" > "+num2str(var_min)+","
				var_test+=1
			endif
			if (var_max<inf&&var_min!=var_max)
				s_constraint+="K"+num2str(v_num_of_waves+1)+" < "+num2str(var_max)+","
				var_test+=1
			endif	
		endif		
	endif
	s_constraint=removeending(s_constraint,",")
	s_constraint+="\""
	//print s_constraint
	if (var_test>0)
		//print s_constraint
		make/D/O M_FitConstraint
		make/D/O W_FitConstraint
		switch (x_shift_done)
			case 1:
				f_constr_txt_to_matrix_MT(s_constraint,(v_num_of_waves+1))
			break
			case 2:
				f_constr_txt_to_matrix_MT(s_constraint,(2*v_num_of_waves+1))
			break
			case 3:
				f_constr_txt_to_matrix_MT(s_constraint,(v_num_of_waves+2))
			break
		
		endswitch
	endif	
	//initial values
	variable var_coef
	switch (x_shift_done)
		case 1: //no x-shift analysis
			make/D/O/N=(v_num_of_waves+1) W_coef_input
			for (index=0;index<(v_num_of_waves+1);index+=1)
				var_coef=start_val[index]
				if (var_coef==0)
					var_coef=1e-6
				endif
				W_coef_input[index]=var_coef
			endfor
		break	
		case 2: //x-shift analysis	multiple waves
			make/D/O/N=(2*v_num_of_waves+1) W_coef_input
			for (index=0;index<(v_num_of_waves+1);index+=1)
				var_coef=start_val[index]
				if (var_coef==0)
					var_coef=1e-6
				endif
				W_coef_input[index]=var_coef
			endfor
			for (index=0;index<(v_num_of_waves);index+=1)	
				var_coef=start_x_val[index]
				W_coef_input[(index+(v_num_of_waves+1))]=var_coef
			endfor		
		break
		case 3:
			make/D/O/N=(v_num_of_waves+2) W_coef_input
			for (index=0;index<(v_num_of_waves+1);index+=1)
				var_coef=start_val[index]
				if (var_coef==0)
					var_coef=1e-6
				endif
				W_coef_input[index]=var_coef
			endfor		
			var_coef=start_x_val[0]
			W_coef_input[(v_num_of_waves+1)]=var_coef
		break
	endswitch
	//correction of start and end values for x-shift analysis
	variable v_start_korr, v_end_korr
	variable v_shift_max_pos,v_shift_max_neg
	switch (x_shift_done)
		case 1:
		break
		case 2: //idea: if the referece is shiftet to negative values, datapoint at the highest values of the wave-to-fit are missing -> fit must ignor these datapoints
			for (index=0;index<(v_num_of_waves);index+=1)
				if (((start_x_val[index])+abs(var_x_wave[index]))>v_shift_max_pos)
					v_shift_max_pos=((start_x_val[index])+abs(var_x_wave[index]))
				endif
				if (((start_x_val[index]-abs(var_x_wave[index])))<v_shift_max_neg)
					v_shift_max_neg=(start_x_val[index]-abs(var_x_wave[index]))
				endif	
				//print (start_x_val[index])+abs(var_x_wave[index])
				//print start_x_val[index]-abs(var_x_wave[index])
			endfor
			v_start_korr=dimoffset(w_sample,2)+v_shift_max_pos+dimdelta(w_sample,2)
			v_end_korr=(dimoffset(w_sample,2)+dimdelta(w_sample,2)*(dimsize(w_sample,2)-1))+v_shift_max_neg-dimdelta(w_sample,2)
			if (v_start<	v_start_korr)
				v_start=v_start_korr
			endif
			if (v_end>	v_end_korr)
				v_end=v_end_korr
			endif		
		break
		case 3:
			v_shift_max_pos=start_x_val[0]+abs(var_x_wave[0])
			v_shift_max_neg=start_x_val[0]-abs(var_x_wave[0])
			//print v_shift_max_pos,v_shift_max_neg
			v_start_korr=dimoffset(w_sample,2)+v_shift_max_pos+dimdelta(w_sample,2)
			v_end_korr=(dimoffset(w_sample,2)+dimdelta(w_sample,2)*(dimsize(w_sample,2)-1))+v_shift_max_neg-dimdelta(w_sample,2)
			//print v_start_korr,v_end_korr
			if (v_start<	v_start_korr)
				v_start=v_start_korr
			endif
			if (v_end>v_end_korr)
				v_end=v_end_korr
			endif					
		break
	endswitch
	variable v_start_energy=v_start
	variable v_end_energy=v_end
	v_start=scaletoindex(w_mask, v_start, 0)
	v_end=scaletoindex(w_mask, v_end, 0)
	//multi-threaded work
	variable fitted_waves
	switch (x_shift_done)
		case 1:
			fitted_waves=v_num_of_waves+1
		break
		case 2:
			fitted_waves=2*v_num_of_waves+1
		break
		case 3:
			fitted_waves=v_num_of_waves+2
		break	
	endswitch
	//print "fitted waves "+num2str(fitted_waves)
	Variable threadGroupID= ThreadGroupCreate(v_threads)
	Variable threadGroupStatus,dummy
	//print v_rows_max
	//print v_lcm_type
	//print "var_test"
	//print var_test
	//print "v_lcm_type"
	//print v_lcm_type
	//print v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim
	if (var_test>0) //only 0 is without contraints
			switch (v_lcm_type)
				case 0:	//point-by-point, no x-shift
					for(t_index=0; t_index<v_threads; t_index+=1)	
						ThreadStart threadGroupID,t_index,f_lc_funcfit_pbp_noxshift_Constraint_MT(w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim)
					endfor	
					do
						threadGroupStatus = ThreadGroupWait(threadGroupID,100)
						while( threadGroupStatus != 0 )
						dummy= ThreadGroupRelease(threadGroupID)
				break
				case 1:	//all-at-ones, no x-shift
					for(t_index=0; t_index<v_threads; t_index+=1)	
						ThreadStart threadGroupID,t_index,f_lc_funcfit_aao_noxshift_Constraint_MT(w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim)
					endfor			
					do
						threadGroupStatus = ThreadGroupWait(threadGroupID,100)
						while( threadGroupStatus != 0 )
						dummy= ThreadGroupRelease(threadGroupID)
				break
				case 2:	//point-by-point, x-shift
					for(t_index=0; t_index<v_threads; t_index+=1)	
						ThreadStart threadGroupID,t_index,f_lc_funcfit_pbp_xshift_Constraint_MT(w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim)
					endfor			
					do
						threadGroupStatus = ThreadGroupWait(threadGroupID,100)
						while( threadGroupStatus != 0 )
						dummy= ThreadGroupRelease(threadGroupID)
				break
				case 3:	//all-at-ones, no x-shift
					for(t_index=0; t_index<v_threads; t_index+=1)	
						ThreadStart threadGroupID,t_index,f_lc_funcfit_aao_xshift_Constraint_MT(w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim)
					endfor			
					do
						threadGroupStatus = ThreadGroupWait(threadGroupID,100)
						while( threadGroupStatus != 0 )
						dummy= ThreadGroupRelease(threadGroupID)
				break
				case 4:
					for(t_index=0; t_index<v_threads; t_index+=1)	
						ThreadStart threadGroupID,t_index,f_lc_funcfit_pbp_single_xshift_Constraint_MT(w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim)
					endfor			
					do
						threadGroupStatus = ThreadGroupWait(threadGroupID,100)
						while( threadGroupStatus != 0 )
						dummy= ThreadGroupRelease(threadGroupID)
				break
				case 5:
					for(t_index=0; t_index<v_threads; t_index+=1)	
						ThreadStart threadGroupID,t_index,f_lc_funcfit_aao_single_xshift_Constraint_MT(w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim)
					endfor			
					do
						threadGroupStatus = ThreadGroupWait(threadGroupID,100)
						while( threadGroupStatus != 0 )
						dummy= ThreadGroupRelease(threadGroupID)	
				break				
			endswitch
		else
			switch (v_lcm_type)
				case 0:	//point-by-point, no x-shift
					for(t_index=0; t_index<v_threads; t_index+=1)	
						ThreadStart threadGroupID,t_index,f_lc_funcfit_pbp_noxshift_noConstraint_MT(w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim)
					endfor	
					do
						threadGroupStatus = ThreadGroupWait(threadGroupID,100)
						while( threadGroupStatus != 0 )
						dummy= ThreadGroupRelease(threadGroupID)
				break
				case 1:	//all-at-ones, no x-shift
					for(t_index=0; t_index<v_threads; t_index+=1)	
						ThreadStart threadGroupID,t_index,f_lc_funcfit_aao_noxshift_noConstraint_MT(w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim)
					endfor			
					do
						threadGroupStatus = ThreadGroupWait(threadGroupID,100)
						while( threadGroupStatus != 0 )
						dummy= ThreadGroupRelease(threadGroupID)
				break
				case 2:	//point-by-point, x-shift, not really existent
					for(t_index=0; t_index<v_threads; t_index+=1)	
						ThreadStart threadGroupID,t_index,f_lc_funcfit_pyp_xshift_noConstraint_MT(w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim)
					endfor			
					do
						threadGroupStatus = ThreadGroupWait(threadGroupID,100)
						while( threadGroupStatus != 0 )
						dummy= ThreadGroupRelease(threadGroupID)
				break
				case 3:	//all-at-ones, no x-shift, not really existent
					for(t_index=0; t_index<v_threads; t_index+=1)	
						ThreadStart threadGroupID,t_index,f_lc_funcfit_aao_xshift_noConstraint_MT(w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim)
					endfor			
					do
						threadGroupStatus = ThreadGroupWait(threadGroupID,100)
						while( threadGroupStatus != 0 )
						dummy= ThreadGroupRelease(threadGroupID)
				break
				case 4:
					for(t_index=0; t_index<v_threads; t_index+=1)	
						ThreadStart threadGroupID,t_index,f_lc_funcfit_pyp_single_xshift_noConstraint_MT(w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim)
					endfor			
					do
						threadGroupStatus = ThreadGroupWait(threadGroupID,100)
						while( threadGroupStatus != 0 )
						dummy= ThreadGroupRelease(threadGroupID)
				break
				case 5:
					for(t_index=0; t_index<v_threads; t_index+=1)	
						ThreadStart threadGroupID,t_index,f_lc_funcfit_aao_single_xshift_noConstraint_MT(w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,fitted_waves,v_LCM_offset,v_LCM_delta, v_fitTol, v_lcm_type, s_hold_LCM, v_inputwave_dim)
					endfor			
					do
						threadGroupStatus = ThreadGroupWait(threadGroupID,100)
						while( threadGroupStatus != 0 )
						dummy= ThreadGroupRelease(threadGroupID)	
				break
			endswitch	
		endif
	//scale
	setscale/P x,dimoffset(w_sample,0),dimdelta(w_sample,0),w_export
	setscale/P x,dimoffset(w_sample,0),dimdelta(w_sample,0),w_stddev
	setscale/P x,dimoffset(w_sample,0),dimdelta(w_sample,0),w_fits
	setscale/P y,dimoffset(w_sample,1),dimdelta(w_sample,1),w_export
	setscale/P y,dimoffset(w_sample,1),dimdelta(w_sample,1),w_stddev
	setscale/P y,dimoffset(w_sample,1),dimdelta(w_sample,1),w_fits
	setscale/P z,dimoffset(w_sample,2),dimdelta(w_sample,2),w_fits		
	//exportfolder
	string s_nameofwave
	string s_datafolder_sample=getwavesDataFolder(w_sample_input,1)
	string exportfolder2 = s_datafolder_sample+"LCMset_"+Secs2Date(DateTime,-2,"_")+"_"									
	string testfolder2
	variable i_folder2
	for (i_folder2=1;i_folder2<inf;i_folder2++)
		testfolder2=exportfolder2+num2str(i_folder2)
		if (datafolderexists(testfolder2)==0)
			newdataFolder $testfolder2
			i_folder2=inf
		endif
	endfor
	setdataFolder $testfolder2
	//generate wave information
	variable v_information_rows
	string s_ref_waves_names=""
	
	if ((v_num_of_waves+2)>10)
		v_information_rows=v_num_of_waves+2
	else
		v_information_rows=10
	endif
	
	switch (x_shift_done)
		case 1:
			make/D/T/O/N=(v_information_rows,5) w_information
		break
		case 2:
			make/D/T/O/N=(v_information_rows,7) w_information
		break
		case 3:
			make/D/T/O/N=(v_information_rows,7) w_information
		break
	endswitch
	w_information[0][0]=Secs2Date(DateTime,2)
	w_information[1][0]=Secs2Time(DateTime,2)
	w_information[2][0]=GetWavesDataFolder(w_sample_input,2)
	if (StringMatch(note(w_sample_input),"")==0)
		w_information[3][0]="Note: "+note(w_sample_input)
	endif
	w_information[6][0]="FitTol: "+num2str(v_fitTol)
	w_information[4][0]="Points: "+num2str(v_start)+" to "+num2str(v_end)
	w_information[5][0]="Scale: "+num2str(v_start_energy)+" to "+num2str(v_end_energy)
	w_information[7][0]="\H="+s_hold_LCM
	w_information[8][0]=s_fit_type
	w_information[9][0]="LCM version "+f_LCM_version()
	w_information[0][1]="References"
	w_information[0][2]="Initial Guess"
	w_information[0][3]="Min"
	w_information[0][4]="Max"
	for (index=0;index<(v_num_of_waves);index++)
		w_information[index+1][1]=GetWavesDataFolder(reference_waves[index],2)
		s_ref_waves_names+=nameofwave(reference_waves[index])+";"
		w_information[index+1][2]=num2str(start_val[index])
		w_information[index+1][3]=num2str(min_wave[index])
		w_information[index+1][4]=num2str(max_wave[index])
	endfor
	w_information[v_num_of_waves+1][1]="absolute Offset"
	s_ref_waves_names+="offset;"
	w_information[v_num_of_waves+1][2]=num2str(start_val[v_num_of_waves])
	w_information[v_num_of_waves+1][3]=num2str(min_wave[v_num_of_waves])
	w_information[v_num_of_waves+1][4]=num2str(max_wave[v_num_of_waves])	
	for (index=0;index<(v_num_of_waves);index++)
		if ((index==0)&&(x_shift_done==3))
			s_ref_waves_names+="Single_xshift;"
		else
			s_ref_waves_names+=nameofwave(reference_waves[index])+"_xshift;"
		endif
	endfor
	//print s_ref_waves_names
	switch (x_shift_done)
	 	case 2:
	 		w_information[0][5]="Initial x-shift"
	 		w_information[0][6]="±"
	 		w_information[8][0]+=" + Multiple x-shift analysis"
			for (index=0;index<v_num_of_waves;index++)
				w_information[index+1][5]=num2str(start_x_val[index])
				w_information[index+1][6]=num2str(var_x_wave[index])
			endfor		 	
	 	break
	 	case 3:
	 		w_information[0][5]="Initial x-shift"
	 		w_information[0][6]="±"
	 		w_information[8][0]+=" + Single x-shift analysis"
	 		w_information[1][5]=num2str(start_x_val[0])
	 		w_information[1][6]=num2str(var_x_wave[0])
	 	break
	endswitch
	//export w_coef & w_sigma
	if (v_inputwave_dim==1)
		f_do_extract_1dim_LCM(w_export, w_stddev ,s_ref_waves_names, nameofwave(w_sample_input), "coef", "sigma")
	else
		if (v_rows==1&&v_columns==1)
			f_do_extract_1dim_LCM(w_export, w_stddev ,s_ref_waves_names, nameofwave(w_sample_input), "coef", "sigma")
		else
			if (v_inputwave_dim==4)
				f_do_extract_4dim_LCM(w_export, s_ref_waves_names, nameofwave(w_sample_input), "coef", v_layers_old)
				f_do_extract_4dim_LCM(w_stddev, s_ref_waves_names, nameofwave(w_sample_input), "sigma", v_layers_old)
			else
				f_do_extract_LCM(w_export, s_ref_waves_names, nameofwave(w_sample_input), "coef")
				f_do_extract_LCM(w_stddev, s_ref_waves_names, nameofwave(w_sample_input), "sigma")
			endif			
		endif
	endif	
	//rename and transfer to subfolder
	string exportname=testfolder2+":"+nameofwave(w_sample_input)+"_coef"
	string exportname2=testfolder2+":"+nameofwave(w_sample_input)+"_sigma"
	string exportname3=testfolder2+":"+"fit_"+nameofwave(w_sample_input)
	string exportname_ref=testfolder2+":"+nameofwave(w_sample_input)+"_references"
	string exportname_mask=testfolder2+":"+nameofwave(w_sample_input)+"_mask"	
	string exportname_info=testfolder2+":"+nameofwave(w_sample_input)+"_information"
	movewave w_export,$exportname
	movewave w_stddev,$exportname2
	movewave w_references_MT,$exportname_ref
	movewave w_mask,$exportname_mask
	movewave w_information,$exportname_info
	setdataFolder $testfolder2	
	killdataFolder/Z root:Packages:LCM:
	f_do_LCM_fit_wave($exportname,$exportname2,$exportname_ref,$exportname_info,v_layers_old,w_sample_input)
	//sound
	if (f_LCM_beep()==1)
		beep
		sleep 00:00:01
		beep
		sleep 00:00:01
		beep
	endif
	print "LCM done."
end

function f_do_LCM_fit_wave(wave w_coef, wave w_sigma, wave w_references, wave/T w_infortmation, variable v_layers_old, wave w_sample_input)
	setdataFolder GetWavesDataFolderDFR(w_coef)
	variable v_waves_max=dimsize(w_references,0)
	variable v_layers=dimsize(w_references,1)
	variable v_coef_max=dimsize(w_coef,2)
	variable v_rows=dimsize(w_coef,0)
	variable v_columns=dimsize(w_coef,1)
	variable v_fit_type=0
	if (v_coef_max==(v_waves_max+1))
		v_fit_type=1
	endif
	if (v_coef_max==(v_waves_max+2))
		v_fit_type=2
	endif
	if (v_coef_max==(2*v_waves_max+1))
		v_fit_type=3
	endif
	variable r_index, c_index, l_index, w_index, v_shift, v_shift_max
	make/FREE/D/O/N=(v_rows,v_columns,v_layers) w_fit=0
	switch (v_fit_type)
		case 1:
			for (r_index=0;r_index<v_rows;r_index++)
				for (c_index=0;c_index<v_columns;c_index++)
					for (w_index=0;w_index<v_waves_max;w_index++)
						for (l_index=0;l_index<v_layers;l_index++)
							w_fit[r_index][c_index][l_index]+=w_coef[r_index][c_index][w_index]*w_references[w_index][l_index]	
						endfor	
					endfor
					for (l_index=0;l_index<v_layers;l_index++)
						w_fit[r_index][c_index][l_index]+=w_coef[r_index][c_index][v_waves_max]
					endfor
				endfor
			endfor			
		break
		case 2:
			for (r_index=0;r_index<v_rows;r_index++)
				for (c_index=0;c_index<v_columns;c_index++)
					for (w_index=0;w_index<v_waves_max;w_index++)
						v_shift=w_coef[r_index][c_index][v_waves_max+1]
						v_shift_max=1.1*abs(str2num(w_infortmation[1][6]))
						if ((v_shift<(-v_shift_max))||(v_shift>v_shift_max))
							for (l_index=0;l_index<v_layers;l_index++)
								w_fit[r_index][c_index][l_index]=NaN
							endfor
						else
							make/FREE/D/O/N=(v_layers) w_temp
							matrixOP/FREE/O w_temp=row(w_references,w_index)
							matrixtranspose w_temp
							setscale/P x,dimoffset(w_references,1),dimdelta(w_references,1),w_temp
							make/FREE/D/O/N=(v_layers) w_temp_scale=dimoffset(w_references,1)+x*dimdelta(w_references,1)
							make/FREE/D/O/N=(v_layers) w_temp_scale_shifted=dimoffset(w_references,1)+x*dimdelta(w_references,1)+v_shift
							for (l_index=0;l_index<v_layers;l_index++)
								w_fit[r_index][c_index][l_index]+=w_coef[r_index][c_index][w_index]*interp((w_temp_scale_shifted[l_index]),w_temp_scale,w_temp)
							endfor	
						endif
					endfor
					for (l_index=0;l_index<v_layers;l_index++)
						w_fit[r_index][c_index][l_index]+=w_coef[r_index][c_index][v_waves_max]
					endfor
				endfor
			endfor			
		break
		case 3:
			for (r_index=0;r_index<v_rows;r_index++)
				for (c_index=0;c_index<v_columns;c_index++)
					for (w_index=0;w_index<v_waves_max;w_index++)
						v_shift=w_coef[r_index][c_index][v_waves_max+1+w_index]
						v_shift_max=1.1*abs(str2num(w_infortmation[1+w_index][6]))
						if ((v_shift<(-v_shift_max))||(v_shift>v_shift_max))
							for (l_index=0;l_index<v_layers;l_index++)
								w_fit[r_index][c_index][l_index]=NaN
							endfor
						else
							make/FREE/D/O/N=(v_layers) w_temp
							matrixOP/FREE/O w_temp=row(w_references,w_index)
							matrixtranspose w_temp
							setscale/P x,dimoffset(w_references,1),dimdelta(w_references,1),w_temp
							make/FREE/D/O/N=(v_layers) w_temp_scale=dimoffset(w_references,1)+x*dimdelta(w_references,1)
							make/FREE/D/O/N=(v_layers) w_temp_scale_shifted=dimoffset(w_references,1)+x*dimdelta(w_references,1)+v_shift
							for (l_index=0;l_index<v_layers;l_index++)
								w_fit[r_index][c_index][l_index]+=w_coef[r_index][c_index][w_index]*interp((w_temp_scale_shifted[l_index]),w_temp_scale,w_temp)
							endfor	
						endif
					endfor
					for (l_index=0;l_index<v_layers;l_index++)
						w_fit[r_index][c_index][l_index]+=w_coef[r_index][c_index][v_waves_max]
					endfor
				endfor
			endfor			
		break
	endswitch
	if ((dimsize(w_fit,0)<2)&&(dimsize(w_fit,1)<2))
		duplicate/FREE w_fit w_fit_3dim
		MatrixOP/FREE/O w_fit=beam(w_fit_3dim,0,0)
		setscale/P x,dimoffset(w_references,1),dimdelta(w_references,1),w_fit
	else	
		setscale/P x,dimoffset(w_coef,0),dimdelta(w_coef,0),w_fit
		setscale/P y,dimoffset(w_coef,1),dimdelta(w_coef,1),w_fit
		setscale/P z,dimoffset(w_references,1),dimdelta(w_references,1),w_fit
	endif
	if (numtype(v_layers_old)==0)
		variable v_chunks=dimsize(w_sample_input,3)
		variable v_columns_new=dimsize(w_sample_input,1)
		v_layers=dimsize(w_sample_input,2)
		v_rows=dimsize(w_sample_input,0)
		make/FREE/N=(v_rows,v_columns_new,v_layers,v_chunks)/O w_4dim
		variable ch_index
		for (r_index=0;r_index<v_rows;r_index++)
			for (c_index=0;c_index<v_columns_new;c_index++)
				for (l_index=0;l_index<v_layers;l_index++)
					for (ch_index=0;ch_index<v_chunks;ch_index++)
						w_4dim[r_index][c_index][l_index][ch_index]=w_fit[r_index][(c_index*v_layers)+l_index][ch_index]
					endfor
				endfor
			endfor	
		endfor	
		setscale/P x,dimoffset(w_sample_input,0),dimdelta(w_sample_input,0),w_4dim	
		setscale/P y,dimoffset(w_sample_input,1),dimdelta(w_sample_input,1),w_4dim	
		setscale/P z,dimoffset(w_sample_input,2),dimdelta(w_sample_input,2),w_4dim	
		setscale/P d,dimoffset(w_sample_input,3),dimdelta(w_sample_input,3),w_4dim	
		duplicate/FREE/O w_4dim w_fit
	endif
	string s_export="fit_"+removeEnding(nameofWave(w_coef),"_coef")
	duplicate/O w_fit $s_export
	killwaves w_coef, w_sigma, w_references //user don't need this waves
end

function f_do_extract_LCM(wave w_matrix ,string s_ref_waves_names, string s_exportname1, string s_exportname2)
	variable v_rows=dimsize(w_matrix,0)
	variable v_columns=dimsize(w_matrix,1)
	variable v_layers=dimsize(w_matrix,2)	
	variable r_index,l_index,c_index, index
	string s_export
	string s_names_ref_waves
	wave wave_ref	
	for (l_index=0;l_index<v_layers;l_index++)
		//s_export=s_exportname1+"_"+StringFromList(l_index,s_ref_waves_names,";")+"_"+s_exportname2  //information of to-be-fitted-wave saved in name
		s_export="LCM_"+StringFromList(l_index,s_ref_waves_names,";")+"_"+s_exportname2  //information of to-be-fitted-wave saved in note
		s_export=replaceString(" ",s_export,"_")
		make/D/N=(v_rows,v_columns) exportwave
		for (r_index=0; r_index<v_rows; r_index +=1)									//change inf to NaN, inf would be automaticly the maximum
			for (c_index=0; c_index<v_columns; c_index +=1)	
				exportwave[r_index][c_index]=w_matrix[r_index][c_index][l_index]
			endfor
		endfor	
		setscale/P x, DimOffset(w_matrix,0), DimDelta(w_matrix,0),exportwave
		setscale/P y, DimOffset(w_matrix,1), DimDelta(w_matrix,1),exportwave	
		note exportwave "LCM of "+s_exportname1 //information of to-be-fitted-wave saved in note
		make/D/O $s_export
		killwaves $s_export
		rename exportwave $s_export	
	endfor
end

function f_do_extract_4dim_LCM(wave w_matrix ,string s_ref_waves_names, string s_exportname1, string s_exportname2, variable v_old_layers)
	variable v_rows=dimsize(w_matrix,0)
	variable v_columns=dimsize(w_matrix,1)
	variable v_layers=dimsize(w_matrix,2)	
	variable v_columns_new=v_columns/v_old_layers
	variable r_index,l_index,c_index,component_index,index
	string s_export
	string s_names_ref_waves
	wave wave_ref	
	for (component_index=0;component_index<v_layers;component_index++)
		//s_export=s_exportname1+"_"+StringFromList(l_index,s_ref_waves_names,";")+"_"+s_exportname2  //information of to-be-fitted-wave saved in name
		s_export="LCM_"+StringFromList(component_index,s_ref_waves_names,";")+"_"+s_exportname2  //information of to-be-fitted-wave saved in note
		s_export=replaceString(" ",s_export,"_")
		make/D/N=(v_rows,v_columns_new,v_old_layers) exportwave
		for (r_index=0; r_index<v_rows; r_index +=1)									//change inf to NaN, inf would be automaticly the maximum
			for (c_index=0; c_index<v_columns_new; c_index +=1)	
				for (l_index=0; l_index<v_old_layers; l_index++)
					exportwave[r_index][c_index][l_index]=w_matrix[r_index][(c_index*v_old_layers)+l_index][component_index]
				endfor
			endfor
		endfor	
		setscale/P x, DimOffset(w_matrix,0), DimDelta(w_matrix,0),exportwave
		setscale/P y, DimOffset(w_matrix,1), DimDelta(w_matrix,1),exportwave	
		note exportwave "LCM of "+s_exportname1 //information of to-be-fitted-wave saved in note
		make/D/O $s_export
		killwaves $s_export
		rename exportwave $s_export	
	endfor
end

function f_do_extract_1dim_LCM(wave w_matrix1, wave w_matrix2 ,string s_ref_waves_names, string s_exportname0, string s_exportname1, string s_exportname2)
	variable v_rows=dimsize(w_matrix1,0)
	variable v_columns=dimsize(w_matrix1,1)
	variable v_layers=dimsize(w_matrix1,2)	
	variable r_index,l_index,c_index, index
	string s_export1, s_export2
	string s_names_ref_waves
	wave wave_ref	
	print "LCM fit for "+s_exportname0+":"
	for (l_index=0;l_index<v_layers;l_index++)
		//s_export1=s_exportname0+"_"+StringFromList(l_index,s_ref_waves_names,";")+"_"+s_exportname1 //information of to-be-fitted-wave saved in name
		//s_export2=s_exportname0+"_"+StringFromList(l_index,s_ref_waves_names,";")+"_"+s_exportname2 //information of to-be-fitted-wave saved in name
		s_export1="LCM_"+StringFromList(l_index,s_ref_waves_names,";")+"_"+s_exportname1  //information of to-be-fitted-wave unsaved
		s_export2="LCM_"+StringFromList(l_index,s_ref_waves_names,";")+"_"+s_exportname2	 //information of to-be-fitted-wave unsaved	
		s_export1=replaceString(" ",s_export1,"_")
		s_export2=replaceString(" ",s_export2,"_")
		variable/G v_export1=w_matrix1[0][0][l_index]
		variable/G v_export2=w_matrix2[0][0][l_index]
		print StringFromList(l_index,s_ref_waves_names,";")+": "+num2str(v_export1)+" ± "+num2str(v_export2)
		rename v_export1 $s_export1
		rename v_export2 $s_export2
	endfor
end



//part of the procedure that is done by each thread
threadsafe function f_lc_funcfit_pbp_noxshift_Constraint_MT(w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, s_hold_LCM, v_inputwave_dim)
 	wave w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT, w_mask
 	variable v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, v_inputwave_dim
 	string s_hold_LCM
 	//points for Thread
 	variable v_rowstart=t_index*v_row_per_thread
	variable v_rowend=(t_index+1)*v_row_per_thread
	if (v_rowend>v_rows_max)
		v_rowend=v_rows_max
	endif
	if (v_rowstart==v_rowend)
		v_rowstart-=1
	endif
	if (v_rowstart<0)
		v_rowstart=0
	endif
	//basic waves	
 	sleep/S t_index+1
 	make/O w_sample_duplicate
 	killwaves w_sample_duplicate
 	duplicate w_sample w_sample_duplicate
 	make/O root:reference
 	killwaves root:reference
 	duplicate w_references_MT root:reference
 	variable row_index,column_index, v_fit_index, v_num_index, index, l_index
 	make/D/O/N=(v_num_of_waves) W_coef 	
 	W_coef=W_coef_input
 	make/D/O/N=(v_num_of_waves) W_sigma
	make/D/N=(v_layers_max)/O tempwave
	make/D/N=(v_layers_max)/O tempwave2
	variable v_real_waves=dimsize(w_references_MT,0)
	make/WAVE/O/N=(v_real_waves) listwave
	string s_name_ref
	//unpack references
	for (index=0; index<(v_real_waves); index+=1)
		make/O/N=(v_layers_max) wave0
		for (l_index=0;l_index<v_layers_max;l_index+=1)
		 wave0[l_index]=w_references_MT[index][l_index]
		endfor
		setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),wave0
		rename wave0 $("reference"+num2str(index))
		listwave[index]=$("root:reference"+num2str(index))
	endfor
	variable/G V_FitTol = V_FitTol_input
	//the fit itself
	for (row_index=v_rowstart;row_index<v_rowend;row_index+=1)
		for (column_index=0;column_index<v_columns_max;column_index+=1)
 				make/D/O/N=(v_num_of_waves) W_coef 	
 				killwaves W_coef
 				duplicate W_coef_input W_coef
 				make/D/O/N=(v_num_of_waves) W_sigma=NaN
				make/D/N=(v_layers_max)/O tempwave=0
				make/D/N=(v_layers_max)/O tempwave2=0
				MatrixOP/O/Z tempwave=beam(w_sample_duplicate,row_index,column_index)
				setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),tempwave,tempwave2
				FuncFit/Q/H=s_hold_LCM f_fitFunc_lc_MT W_coef tempwave[v_start,v_end] /C={M_FitConstraint,W_FitConstraint} /M=w_mask /NWOK
				//save coefficients
				for (v_num_index=0;v_num_index<(v_num_of_waves);v_num_index+=1)
					w_export[row_index][column_index][v_num_index]=W_coef[v_num_index]
					w_stddev[row_index][column_index][v_num_index]=W_sigma[v_num_index]
				endfor
			endfor
		endfor
end

//part of the procedure that is done by each thread
threadsafe function f_lc_funcfit_aao_noxshift_Constraint_MT(w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, s_hold_LCM, v_inputwave_dim)
 	wave w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT, w_mask
 	variable v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, v_inputwave_dim
 	string s_hold_LCM
 	//points for Thread
 	variable v_rowstart=t_index*v_row_per_thread
	variable v_rowend=(t_index+1)*v_row_per_thread
	if (v_rowend>v_rows_max)
		v_rowend=v_rows_max
	endif
	if (v_rowstart>v_rows_max)
		v_rowstart=v_rows_max
	endif
	if (v_rowstart>v_rowend)
		v_rowstart=v_rowend-1
	endif		
	if (v_rowstart==v_rowend)
		v_rowstart-=1
	endif
	if (v_rowstart<0)
		v_rowstart=0
	endif
 	sleep/S t_index+1
 	make/O w_sample_duplicate
 	killwaves w_sample_duplicate
 	duplicate w_sample w_sample_duplicate
 	make/O root:reference
 	killwaves root:reference
 	duplicate w_references_MT root:reference
 	variable row_index,column_index, v_fit_index, v_num_index, index, l_index
 	make/D/O/N=(v_num_of_waves) W_coef 	
 	W_coef=W_coef_input
 	make/D/O/N=(v_num_of_waves) W_sigma
	make/D/N=(v_layers_max)/O tempwave
	make/D/N=(v_layers_max)/O tempwave2
	variable v_real_waves=dimsize(w_references_MT,0)
	make/WAVE/O/N=(v_real_waves) listwave
	string s_name_ref
	//unpack references
	for (index=0; index<(v_real_waves); index+=1)
		make/O/N=(v_layers_max) wave0
		for (l_index=0;l_index<v_layers_max;l_index+=1)
		 wave0[l_index]=w_references_MT[index][l_index]
		endfor
		setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),wave0
		rename wave0 $("reference"+num2str(index))
		listwave[index]=$("root:reference"+num2str(index))
	endfor
	variable/G V_FitTol = V_FitTol_input   
	//the fit itself
	for (row_index=v_rowstart;row_index<v_rowend;row_index+=1)
		for (column_index=0;column_index<v_columns_max;column_index+=1)
 				make/D/O/N=(v_num_of_waves) W_coef 	
 				killwaves W_coef
 				duplicate W_coef_input W_coef
 				make/D/O/N=(v_num_of_waves) W_sigma=NaN
				make/D/N=(v_layers_max)/O tempwave=0
				make/D/N=(v_layers_max)/O tempwave2=0
				MatrixOP/O/Z tempwave=beam(w_sample_duplicate,row_index,column_index)
				setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),tempwave,tempwave2
				FuncFit/Q/H=s_hold_LCM f_fitFunc_lc_MT_AllAtOnce W_coef tempwave[v_start,v_end] /C={M_FitConstraint,W_FitConstraint} /M=w_mask /NWOK
				//save coefficients
				for (v_num_index=0;v_num_index<(v_num_of_waves);v_num_index+=1)
					w_export[row_index][column_index][v_num_index]=W_coef[v_num_index]
					w_stddev[row_index][column_index][v_num_index]=W_sigma[v_num_index]
				endfor
			endfor
		endfor
end

//part of the procedure that is done by each thread
threadsafe function f_lc_funcfit_pbp_xshift_Constraint_MT(w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, s_hold_LCM, v_inputwave_dim)
 	wave w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT, w_mask
 	variable v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, v_inputwave_dim
 	string s_hold_LCM
 	//points for Thread
 	variable v_rowstart=t_index*v_row_per_thread
	variable v_rowend=(t_index+1)*v_row_per_thread
	if (v_rowend>v_rows_max)
		v_rowend=v_rows_max
	endif
	if (v_rowstart>v_rows_max)
		v_rowstart=v_rows_max
	endif
	if (v_rowstart>v_rowend)
		v_rowstart=v_rowend-1
	endif		
	if (v_rowstart==v_rowend)
		v_rowstart-=1
	endif
	if (v_rowstart<0)
		v_rowstart=0
	endif
	//basic waves	
 	sleep/S t_index+1
 	make/O w_sample_duplicate
 	killwaves w_sample_duplicate
 	duplicate w_sample w_sample_duplicate
 	make/O root:reference
 	killwaves root:reference
 	duplicate w_references_MT root:reference
 	variable row_index,column_index, v_fit_index, v_num_index, index, l_index
 	make/D/O/N=(v_num_of_waves) W_coef 	
 	W_coef=W_coef_input
 	make/D/O/N=(v_num_of_waves) W_sigma
	make/D/N=(v_layers_max)/O tempwave
	make/D/N=(v_layers_max)/O tempwave2
	variable v_real_waves=dimsize(w_references_MT,0)
	make/WAVE/O/N=(v_real_waves) listwave
	string s_name_ref
	//unpack references
	for (index=0; index<(v_real_waves); index+=1)
		make/O/N=(v_layers_max) wave0
		for (l_index=0;l_index<v_layers_max;l_index+=1)
		 wave0[l_index]=w_references_MT[index][l_index]
		endfor
		setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),wave0
		rename wave0 $("reference"+num2str(index))
		listwave[index]=$("root:reference"+num2str(index))
	endfor
	variable/G V_FitTol = V_FitTol_input
	//the fit itself
	for (row_index=v_rowstart;row_index<v_rowend;row_index+=1)
		for (column_index=0;column_index<v_columns_max;column_index+=1)
 			make/D/O/N=(v_num_of_waves) W_coef 	
 			killwaves W_coef
 			duplicate W_coef_input W_coef
 			make/D/O/N=(v_num_of_waves) W_sigma=NaN
			make/D/N=(v_layers_max)/O tempwave=0
			make/D/N=(v_layers_max)/O tempwave2=0
			MatrixOP/O/Z tempwave=beam(w_sample_duplicate,row_index,column_index)
			setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),tempwave,tempwave2
			FuncFit/Q/H=s_hold_LCM  f_fitFunc_lc_MT_multishift W_coef tempwave[v_start,v_end] /C={M_FitConstraint,W_FitConstraint} /M=w_mask /NWOK // /H=s_hold_LCM 		
			//save coefficients
			for (v_num_index=0;v_num_index<(v_num_of_waves);v_num_index+=1)
				w_export[row_index][column_index][v_num_index]=W_coef[v_num_index]
				w_stddev[row_index][column_index][v_num_index]=W_sigma[v_num_index]
			endfor
		endfor
	endfor
end


threadsafe function f_lc_funcfit_pbp_single_xshift_Constraint_MT(w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, s_hold_LCM, v_inputwave_dim)
 	wave w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT, w_mask
 	variable v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, v_inputwave_dim
 	string s_hold_LCM
 	//points for Thread
 	variable v_rowstart=t_index*v_row_per_thread
	variable v_rowend=(t_index+1)*v_row_per_thread
	if (v_rowend>v_rows_max)
		v_rowend=v_rows_max
	endif
	if (v_rowstart>v_rows_max)
		v_rowstart=v_rows_max
	endif
	if (v_rowstart>v_rowend)
		v_rowstart=v_rowend-1
	endif		
	if (v_rowstart==v_rowend)
		v_rowstart-=1
	endif
	if (v_rowstart<0)
		v_rowstart=0
	endif
	//basic waves	
 	sleep/S t_index+1
 	make/O w_sample_duplicate
 	killwaves w_sample_duplicate
 	duplicate w_sample w_sample_duplicate
 	make/O root:reference
 	killwaves root:reference
 	duplicate w_references_MT root:reference
 	variable row_index,column_index, v_fit_index, v_num_index, index, l_index
 	make/D/O/N=(v_num_of_waves) W_coef 	
 	W_coef=W_coef_input
 	make/D/O/N=(v_num_of_waves) W_sigma
	make/D/N=(v_layers_max)/O tempwave
	make/D/N=(v_layers_max)/O tempwave2
	variable v_real_waves=dimsize(w_references_MT,0)
	make/WAVE/O/N=(v_real_waves) listwave
	string s_name_ref
	//unpack references
	for (index=0; index<(v_real_waves); index+=1)
		make/O/N=(v_layers_max) wave0
		for (l_index=0;l_index<v_layers_max;l_index+=1)
		 wave0[l_index]=w_references_MT[index][l_index]
		endfor
		setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),wave0
		rename wave0 $("reference"+num2str(index))
		listwave[index]=$("root:reference"+num2str(index))
	endfor
	variable/G V_FitTol = V_FitTol_input
	//the fit itself
	for (row_index=v_rowstart;row_index<v_rowend;row_index+=1)
		for (column_index=0;column_index<v_columns_max;column_index+=1)
 			make/D/O/N=(v_num_of_waves) W_coef 	
 			killwaves W_coef
 			duplicate W_coef_input W_coef
 			make/D/O/N=(v_num_of_waves) W_sigma=NaN
			make/D/N=(v_layers_max)/O tempwave=0
			make/D/N=(v_layers_max)/O tempwave2=0
			MatrixOP/O/Z tempwave=beam(w_sample_duplicate,row_index,column_index)
			setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),tempwave,tempwave2
			FuncFit/Q/H=s_hold_LCM  f_fitFunc_lc_MT_singleshift W_coef tempwave[v_start,v_end] /C={M_FitConstraint,W_FitConstraint} /M=w_mask /NWOK // /H=s_hold_LCM 		
			//save coefficients
			for (v_num_index=0;v_num_index<(v_num_of_waves);v_num_index+=1)
				w_export[row_index][column_index][v_num_index]=W_coef[v_num_index]
				w_stddev[row_index][column_index][v_num_index]=W_sigma[v_num_index]
			endfor
		endfor
	endfor
end


//part of the procedure that is done by each thread
threadsafe function f_lc_funcfit_aao_xshift_Constraint_MT(w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, s_hold_LCM, v_inputwave_dim)
 	wave w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT, w_mask
 	variable v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, v_inputwave_dim
 	string s_hold_LCM
 	//points for Thread
 	variable v_rowstart=t_index*v_row_per_thread
	variable v_rowend=(t_index+1)*v_row_per_thread
	if (v_rowend>v_rows_max)
		v_rowend=v_rows_max
	endif
	if (v_rowstart>v_rows_max)
		v_rowstart=v_rows_max
	endif
	if (v_rowstart>v_rowend)
		v_rowstart=v_rowend-1
	endif		
	if (v_rowstart==v_rowend)
		v_rowstart-=1
	endif
	if (v_rowstart<0)
		v_rowstart=0
	endif
	//basic waves	
 	sleep/S t_index+1
 	make/O w_sample_duplicate
 	killwaves w_sample_duplicate
 	duplicate w_sample w_sample_duplicate
 	make/O root:reference
 	killwaves root:reference
 	duplicate w_references_MT root:reference
 	variable row_index,column_index, v_fit_index, v_num_index, index, l_index
 	make/D/O/N=(v_num_of_waves) W_coef 	
 	W_coef=W_coef_input
 	make/D/O/N=(v_num_of_waves) W_sigma
	make/D/N=(v_layers_max)/O tempwave
	make/D/N=(v_layers_max)/O tempwave2
	variable v_real_waves=dimsize(w_references_MT,0)
	make/WAVE/O/N=(v_real_waves) listwave
	string s_name_ref
	//unpack references
	for (index=0; index<(v_real_waves); index+=1)
		make/O/N=(v_layers_max) wave0
		for (l_index=0;l_index<v_layers_max;l_index+=1)
		 wave0[l_index]=w_references_MT[index][l_index]
		endfor
		setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),wave0
		rename wave0 $("reference"+num2str(index))
		listwave[index]=$("root:reference"+num2str(index))
	endfor
	variable/G V_FitTol = V_FitTol_input
	//the fit itself
	for (row_index=v_rowstart;row_index<v_rowend;row_index+=1)
		for (column_index=0;column_index<v_columns_max;column_index+=1)
 			make/D/O/N=(v_num_of_waves) W_coef 	
 			killwaves W_coef
 			duplicate W_coef_input W_coef
 			make/D/O/N=(v_num_of_waves) W_sigma=NaN
			make/D/N=(v_layers_max)/O tempwave=0
			make/D/N=(v_layers_max)/O tempwave2=0
			MatrixOP/O/Z tempwave=beam(w_sample_duplicate,row_index,column_index)
			setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),tempwave,tempwave2
			FuncFit/Q/H=s_hold_LCM  f_fitFunc_lc_MT_AllAtOnce_multishift W_coef tempwave[v_start,v_end] /C={M_FitConstraint,W_FitConstraint} /M=w_mask /NWOK  // /H=s_hold_LCM 		
			//save coefficients
			for (v_num_index=0;v_num_index<(v_num_of_waves);v_num_index+=1)
				w_export[row_index][column_index][v_num_index]=W_coef[v_num_index]
				w_stddev[row_index][column_index][v_num_index]=W_sigma[v_num_index]
			endfor
		endfor
	endfor
end

threadsafe function f_lc_funcfit_aao_single_xshift_Constraint_MT(w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, s_hold_LCM, v_inputwave_dim)
 	wave w_sample,M_FitConstraint,W_FitConstraint,W_coef_input,w_export,w_stddev,w_fits,w_references_MT, w_mask
 	variable v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, v_inputwave_dim
 	string s_hold_LCM
 	//points for Thread
 	variable v_rowstart=t_index*v_row_per_thread
	variable v_rowend=(t_index+1)*v_row_per_thread
	if (v_rowend>v_rows_max)
		v_rowend=v_rows_max
	endif
	if (v_rowstart>v_rows_max)
		v_rowstart=v_rows_max
	endif
	if (v_rowstart>v_rowend)
		v_rowstart=v_rowend-1
	endif		
	if (v_rowstart==v_rowend)
		v_rowstart-=1
	endif
	if (v_rowstart<0)
		v_rowstart=0
	endif
	//basic waves	
 	sleep/S t_index+1
 	make/O w_sample_duplicate
 	killwaves w_sample_duplicate
 	duplicate w_sample w_sample_duplicate
 	make/O root:reference
 	killwaves root:reference
 	duplicate w_references_MT root:reference
 	variable row_index,column_index, v_fit_index, v_num_index, index, l_index
 	make/D/O/N=(v_num_of_waves) W_coef 	
 	W_coef=W_coef_input
 	make/D/O/N=(v_num_of_waves) W_sigma
	make/D/N=(v_layers_max)/O tempwave
	make/D/N=(v_layers_max)/O tempwave2
	variable v_real_waves=dimsize(w_references_MT,0)
	make/WAVE/O/N=(v_real_waves) listwave
	string s_name_ref
	//unpack references
	for (index=0; index<(v_real_waves); index+=1)
		make/O/N=(v_layers_max) wave0
		for (l_index=0;l_index<v_layers_max;l_index+=1)
		 wave0[l_index]=w_references_MT[index][l_index]
		endfor
		setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),wave0
		rename wave0 $("reference"+num2str(index))
		listwave[index]=$("root:reference"+num2str(index))
	endfor
	variable/G V_FitTol = V_FitTol_input
	//the fit itself
	for (row_index=v_rowstart;row_index<v_rowend;row_index+=1)
		for (column_index=0;column_index<v_columns_max;column_index+=1)
 			make/D/O/N=(v_num_of_waves) W_coef 	
 			killwaves W_coef
 			duplicate W_coef_input W_coef
 			make/D/O/N=(v_num_of_waves) W_sigma=NaN
			make/D/N=(v_layers_max)/O tempwave=0
			make/D/N=(v_layers_max)/O tempwave2=0
			MatrixOP/O/Z tempwave=beam(w_sample_duplicate,row_index,column_index)
			setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),tempwave,tempwave2
			FuncFit/Q/H=s_hold_LCM  f_fitFunc_lc_MT_AllAtOnce_singleshift W_coef tempwave[v_start,v_end] /C={M_FitConstraint,W_FitConstraint} /M=w_mask /NWOK  // /H=s_hold_LCM 		
			//save coefficients
			for (v_num_index=0;v_num_index<(v_num_of_waves);v_num_index+=1)
				w_export[row_index][column_index][v_num_index]=W_coef[v_num_index]
				w_stddev[row_index][column_index][v_num_index]=W_sigma[v_num_index]
			endfor
		endfor
	endfor
end

threadsafe function f_lc_funcfit_pbp_noxshift_noConstraint_MT(w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, s_hold_LCM, v_inputwave_dim)
 	wave w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask
 	variable v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, v_inputwave_dim
 	string s_hold_LCM
 	//print v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input
 	//points for Thread
 	variable v_rowstart=t_index*v_row_per_thread
	variable v_rowend=(t_index+1)*v_row_per_thread
	if (v_rowend>v_rows_max)
		v_rowend=v_rows_max
	endif
	if (v_rowstart>v_rows_max)
		v_rowstart=v_rows_max-1
	endif
	if (v_rowstart>v_rowend)
		v_rowstart=v_rowend-1
	endif		
	if (v_rowstart==v_rowend)
		v_rowstart-=1
	endif
	if (v_rowstart<0)
		v_rowstart=0
	endif
	//basic waves	
 	sleep/S t_index+1
 	make/O w_sample_duplicate
 	killwaves w_sample_duplicate
 	duplicate w_sample w_sample_duplicate
 	make/O root:reference
 	killwaves root:reference
 	duplicate w_references_MT root:reference
 	variable row_index,column_index, v_fit_index, v_num_index, index, l_index
 	make/D/O/N=(v_num_of_waves) W_coef 	
 	W_coef=W_coef_input
 	make/D/O/N=(v_num_of_waves) W_sigma
	make/D/N=(v_layers_max)/O tempwave
	make/D/N=(v_layers_max)/O tempwave2
	variable v_real_waves=dimsize(w_references_MT,0)
	make/WAVE/O/N=(v_real_waves) listwave
	string s_name_ref
	//unpack references
	for (index=0; index<(v_real_waves); index+=1)
		make/O/N=(v_layers_max) wave0
		for (l_index=0;l_index<v_layers_max;l_index+=1)
		 wave0[l_index]=w_references_MT[index][l_index]
		endfor
		setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),wave0
		rename wave0 $("reference"+num2str(index))
		listwave[index]=$("root:reference"+num2str(index))
	endfor
	variable/G V_FitTol = V_FitTol_input
	//the fit itself
	for (row_index=v_rowstart;row_index<v_rowend;row_index+=1)
		for (column_index=0;column_index<v_columns_max;column_index+=1)
 			make/D/O/N=(v_num_of_waves) W_coef 	
 			killwaves W_coef
 			duplicate W_coef_input W_coef
 			make/D/O/N=(v_num_of_waves) W_sigma=NaN
			make/D/N=(v_layers_max)/O tempwave=0
			make/D/N=(v_layers_max)/O tempwave2=0
			MatrixOP/O/Z tempwave=beam(w_sample_duplicate,row_index,column_index)
			setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),tempwave,tempwave2
			FuncFit/Q/H=s_hold_LCM f_fitFunc_lc_MT W_coef tempwave[v_start,v_end] /M=w_mask /NWOK
			//save coefficients
			for (v_num_index=0;v_num_index<(v_num_of_waves);v_num_index+=1)
				w_export[row_index][column_index][v_num_index]=W_coef[v_num_index]
				w_stddev[row_index][column_index][v_num_index]=W_sigma[v_num_index]
			endfor
		endfor
	endfor
end

threadsafe function f_lc_funcfit_aao_noxshift_noConstraint_MT(w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, s_hold_LCM, v_inputwave_dim)
 	wave w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask
 	variable v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, v_inputwave_dim
 	string s_hold_LCM
 	//print v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input
 	//points for Thread
 	variable v_rowstart=t_index*v_row_per_thread
	variable v_rowend=(t_index+1)*v_row_per_thread
	if (v_rowend>v_rows_max)
		v_rowend=v_rows_max
	endif
	if (v_rowstart>v_rows_max)
		v_rowstart=v_rows_max
	endif
	if (v_rowstart>v_rowend)
		v_rowstart=v_rowend-1
	endif		
	if (v_rowstart==v_rowend)
		v_rowstart-=1
	endif
	if (v_rowstart<0)
		v_rowstart=0
	endif
	//print "v_rowstart: "+num2str(v_rowstart)
	//print "v_rowend: "+num2str(v_rowend)
	//basic waves	
 	sleep/S t_index+1
 	make/O w_sample_duplicate
 	killwaves w_sample_duplicate
 	duplicate w_sample w_sample_duplicate
 	make/O root:reference
 	killwaves root:reference
 	duplicate w_references_MT root:reference
 	variable row_index,column_index, v_fit_index, v_num_index, index, l_index
 	make/D/O/N=(v_num_of_waves) W_coef 	
 	W_coef=W_coef_input
 	make/D/O/N=(v_num_of_waves) W_sigma
	make/D/N=(v_layers_max)/O tempwave
	make/D/N=(v_layers_max)/O tempwave2
	variable v_real_waves=dimsize(w_references_MT,0)
	make/WAVE/O/N=(v_real_waves) listwave
	string s_name_ref
	//unpack references
	for (index=0; index<(v_real_waves); index+=1)
		make/O/N=(v_layers_max) wave0
		for (l_index=0;l_index<v_layers_max;l_index+=1)
		 wave0[l_index]=w_references_MT[index][l_index]
		endfor
		setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),wave0
		rename wave0 $("reference"+num2str(index))
		listwave[index]=$("root:reference"+num2str(index))
	endfor
	variable/G V_FitTol = V_FitTol_input
	//the fit itself
	for (row_index=v_rowstart;row_index<v_rowend;row_index+=1)
		//print row_index
		for (column_index=0;column_index<v_columns_max;column_index+=1)
 			make/D/O/N=(v_num_of_waves) W_coef 	
 			killwaves W_coef
 			duplicate W_coef_input W_coef
 			make/D/O/N=(v_num_of_waves) W_sigma=NaN
			make/D/N=(v_layers_max)/O tempwave=0
			make/D/N=(v_layers_max)/O tempwave2=0
			MatrixOP/O/Z tempwave=beam(w_sample_duplicate,row_index,column_index)
			setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),tempwave,tempwave2
			FuncFit/Q/H=s_hold_LCM f_fitFunc_lc_MT_AllAtOnce W_coef tempwave[v_start,v_end] /M=w_mask /NWOK
			//save coefficients
			for (v_num_index=0;v_num_index<(v_num_of_waves);v_num_index+=1)
				w_export[row_index][column_index][v_num_index]=W_coef[v_num_index]
				w_stddev[row_index][column_index][v_num_index]=W_sigma[v_num_index]
			endfor
		endfor
	endfor
end

threadsafe function f_lc_funcfit_pyp_xshift_noConstraint_MT(w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, s_hold_LCM, v_inputwave_dim)
 	wave w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask
 	variable v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, v_inputwave_dim
 	string s_hold_LCM
 	//print v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input
 	//points for Thread
 	variable v_rowstart=t_index*v_row_per_thread
	variable v_rowend=(t_index+1)*v_row_per_thread
	if (v_rowend>v_rows_max)
		v_rowend=v_rows_max
	endif
	if (v_rowstart>v_rows_max)
		v_rowstart=v_rows_max
	endif
	if (v_rowstart>v_rowend)
		v_rowstart=v_rowend-1
	endif		
	if (v_rowstart==v_rowend)
		v_rowstart-=1
	endif
	if (v_rowstart<0)
		v_rowstart=0
	endif
	//basic waves	
 	sleep/S t_index+1
 	make/O w_sample_duplicate
 	killwaves w_sample_duplicate
 	duplicate w_sample w_sample_duplicate
 	make/O root:reference
 	killwaves root:reference
 	duplicate w_references_MT root:reference
 	variable row_index,column_index, v_fit_index, v_num_index, index, l_index
 	make/D/O/N=(v_num_of_waves) W_coef 	
 	W_coef=W_coef_input
 	make/D/O/N=(v_num_of_waves) W_sigma
	make/D/N=(v_layers_max)/O tempwave
	make/D/N=(v_layers_max)/O tempwave2
	variable v_real_waves=dimsize(w_references_MT,0)
	make/WAVE/O/N=(v_real_waves) listwave
	string s_name_ref
	//unpack references
	for (index=0; index<(v_real_waves); index+=1)
		make/O/N=(v_layers_max) wave0
		for (l_index=0;l_index<v_layers_max;l_index+=1)
		 wave0[l_index]=w_references_MT[index][l_index]
		endfor
		setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),wave0
		rename wave0 $("reference"+num2str(index))
		listwave[index]=$("root:reference"+num2str(index))
	endfor
	variable/G V_FitTol = V_FitTol_input
	//the fit itself
	for (row_index=v_rowstart;row_index<v_rowend;row_index+=1)
		for (column_index=0;column_index<v_columns_max;column_index+=1)
 			make/D/O/N=(v_num_of_waves) W_coef 	
 			killwaves W_coef
 			duplicate W_coef_input W_coef
 			make/D/O/N=(v_num_of_waves) W_sigma=NaN
			make/D/N=(v_layers_max)/O tempwave=0
			make/D/N=(v_layers_max)/O tempwave2=0
			MatrixOP/O/Z tempwave=beam(w_sample_duplicate,row_index,column_index)
			setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),tempwave,tempwave2
			FuncFit/Q/H=s_hold_LCM   f_fitFunc_lc_MT_multishift W_coef tempwave[v_start,v_end] /M=w_mask /NWOK  // /H=s_hold_LCM 
			//save coefficients
			for (v_num_index=0;v_num_index<(v_num_of_waves);v_num_index+=1)
				w_export[row_index][column_index][v_num_index]=W_coef[v_num_index]
				w_stddev[row_index][column_index][v_num_index]=W_sigma[v_num_index]
			endfor
		endfor
	endfor
end

threadsafe function f_lc_funcfit_pyp_single_xshift_noConstraint_MT(w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, s_hold_LCM, v_inputwave_dim)
 	wave w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask
 	variable v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, v_inputwave_dim
 	string s_hold_LCM
 	//print v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input
 	//points for Thread
 	variable v_rowstart=t_index*v_row_per_thread
	variable v_rowend=(t_index+1)*v_row_per_thread
	if (v_rowend>v_rows_max)
		v_rowend=v_rows_max
	endif
	if (v_rowstart>v_rows_max)
		v_rowstart=v_rows_max
	endif
	if (v_rowstart>v_rowend)
		v_rowstart=v_rowend-1
	endif		
	if (v_rowstart==v_rowend)
		v_rowstart-=1
	endif
	if (v_rowstart<0)
		v_rowstart=0
	endif
	//basic waves	
 	sleep/S t_index+1
 	make/O w_sample_duplicate
 	killwaves w_sample_duplicate
 	duplicate w_sample w_sample_duplicate
 	make/O root:reference
 	killwaves root:reference
 	duplicate w_references_MT root:reference
 	variable row_index,column_index, v_fit_index, v_num_index, index, l_index
 	make/D/O/N=(v_num_of_waves) W_coef 	
 	W_coef=W_coef_input
 	make/D/O/N=(v_num_of_waves) W_sigma
	make/D/N=(v_layers_max)/O tempwave
	make/D/N=(v_layers_max)/O tempwave2
	variable v_real_waves=dimsize(w_references_MT,0)
	make/WAVE/O/N=(v_real_waves) listwave
	string s_name_ref
	//unpack references
	for (index=0; index<(v_real_waves); index+=1)
		make/O/N=(v_layers_max) wave0
		for (l_index=0;l_index<v_layers_max;l_index+=1)
		 wave0[l_index]=w_references_MT[index][l_index]
		endfor
		setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),wave0
		rename wave0 $("reference"+num2str(index))
		listwave[index]=$("root:reference"+num2str(index))
	endfor
	variable/G V_FitTol = V_FitTol_input
	//the fit itself
	for (row_index=v_rowstart;row_index<v_rowend;row_index+=1)
		for (column_index=0;column_index<v_columns_max;column_index+=1)
 			make/D/O/N=(v_num_of_waves) W_coef 	
 			killwaves W_coef
 			duplicate W_coef_input W_coef
 			make/D/O/N=(v_num_of_waves) W_sigma=NaN
			make/D/N=(v_layers_max)/O tempwave=0
			make/D/N=(v_layers_max)/O tempwave2=0
			MatrixOP/O/Z tempwave=beam(w_sample_duplicate,row_index,column_index)
			setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),tempwave,tempwave2
			FuncFit/Q/H=s_hold_LCM   f_fitFunc_lc_MT_singleshift W_coef tempwave[v_start,v_end] /M=w_mask /NWOK  // /H=s_hold_LCM 
			//save coefficients
			for (v_num_index=0;v_num_index<(v_num_of_waves);v_num_index+=1)
				w_export[row_index][column_index][v_num_index]=W_coef[v_num_index]
				w_stddev[row_index][column_index][v_num_index]=W_sigma[v_num_index]
			endfor
		endfor
	endfor
end

threadsafe function f_lc_funcfit_aao_xshift_noConstraint_MT(w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, s_hold_LCM, v_inputwave_dim)
 	wave w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask
 	variable v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, v_inputwave_dim
 	string s_hold_LCM
 	//print v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input
 	//points for Thread
 	variable v_rowstart=t_index*v_row_per_thread
	variable v_rowend=(t_index+1)*v_row_per_thread
	if (v_rowend>v_rows_max)
		v_rowend=v_rows_max
	endif
	if (v_rowstart>v_rows_max)
		v_rowstart=v_rows_max
	endif
	if (v_rowstart>v_rowend)
		v_rowstart=v_rowend-1
	endif		
	if (v_rowstart==v_rowend)
		v_rowstart-=1
	endif
	if (v_rowstart<0)
		v_rowstart=0
	endif
	//basic waves	
 	sleep/S t_index+1
 	make/O w_sample_duplicate
 	killwaves w_sample_duplicate
 	duplicate w_sample w_sample_duplicate
 	make/O root:reference
 	killwaves root:reference
 	duplicate w_references_MT root:reference
 	variable row_index,column_index, v_fit_index, v_num_index, index, l_index
 	make/D/O/N=(v_num_of_waves) W_coef 	
 	W_coef=W_coef_input
 	make/D/O/N=(v_num_of_waves) W_sigma
	make/D/N=(v_layers_max)/O tempwave
	make/D/N=(v_layers_max)/O tempwave2
	variable v_real_waves=dimsize(w_references_MT,0)
	make/WAVE/O/N=(v_real_waves) listwave
	string s_name_ref
	//unpack references
	for (index=0; index<(v_real_waves); index+=1)
		make/O/N=(v_layers_max) wave0
		for (l_index=0;l_index<v_layers_max;l_index+=1)
		 wave0[l_index]=w_references_MT[index][l_index]
		endfor
		setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),wave0
		rename wave0 $("reference"+num2str(index))
		listwave[index]=$("root:reference"+num2str(index))
	endfor
	variable/G V_FitTol = V_FitTol_input
	//the fit itself
	for (row_index=v_rowstart;row_index<v_rowend;row_index+=1)
		for (column_index=0;column_index<v_columns_max;column_index+=1)
 			make/D/O/N=(v_num_of_waves) W_coef 	
 			killwaves W_coef
 			duplicate W_coef_input W_coef
 			make/D/O/N=(v_num_of_waves) W_sigma=NaN
			make/D/N=(v_layers_max)/O tempwave=0
			make/D/N=(v_layers_max)/O tempwave2=0
			MatrixOP/O/Z tempwave=beam(w_sample_duplicate,row_index,column_index)
			setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),tempwave,tempwave2
			FuncFit/Q/H=s_hold_LCM   f_fitFunc_lc_MT_AllAtOnce_multishift W_coef tempwave[v_start,v_end] /M=w_mask /NWOK // /H=s_hold_LCM 
			//save coefficients
			for (v_num_index=0;v_num_index<(v_num_of_waves);v_num_index+=1)
				w_export[row_index][column_index][v_num_index]=W_coef[v_num_index]
				w_stddev[row_index][column_index][v_num_index]=W_sigma[v_num_index]
			endfor

		endfor
	endfor
end

threadsafe function f_lc_funcfit_aao_single_xshift_noConstraint_MT(w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask,v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, s_hold_LCM, v_inputwave_dim)
 	wave w_sample,W_coef_input,w_export,w_stddev,w_fits,w_references_MT,w_mask
 	variable v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input, v_lcm_type, v_inputwave_dim
 	string s_hold_LCM
 	//print v_rows_max,v_columns_max,v_layers_max,v_row_per_thread,t_index,v_start,v_end,v_num_of_waves,v_sampleoffset,v_sampledelta, V_FitTol_input
 	//points for Thread
 	variable v_rowstart=t_index*v_row_per_thread
	variable v_rowend=(t_index+1)*v_row_per_thread
	if (v_rowend>v_rows_max)
		v_rowend=v_rows_max
	endif
	if (v_rowstart>v_rows_max)
		v_rowstart=v_rows_max
	endif
	if (v_rowstart>v_rowend)
		v_rowstart=v_rowend-1
	endif		
	if (v_rowstart==v_rowend)
		v_rowstart-=1
	endif
	if (v_rowstart<0)
		v_rowstart=0
	endif
	//basic waves	
 	sleep/S t_index+1
 	make/O w_sample_duplicate
 	killwaves w_sample_duplicate
 	duplicate w_sample w_sample_duplicate
 	make/O root:reference
 	killwaves root:reference
 	duplicate w_references_MT root:reference
 	variable row_index,column_index, v_fit_index, v_num_index, index, l_index
 	make/D/O/N=(v_num_of_waves) W_coef 	
 	W_coef=W_coef_input
 	make/D/O/N=(v_num_of_waves) W_sigma
	make/D/N=(v_layers_max)/O tempwave
	make/D/N=(v_layers_max)/O tempwave2
	variable v_real_waves=dimsize(w_references_MT,0)
	make/WAVE/O/N=(v_real_waves) listwave
	string s_name_ref
	//unpack references
	for (index=0; index<(v_real_waves); index+=1)
		make/O/N=(v_layers_max) wave0
		for (l_index=0;l_index<v_layers_max;l_index+=1)
		 wave0[l_index]=w_references_MT[index][l_index]
		endfor
		setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),wave0
		rename wave0 $("reference"+num2str(index))
		listwave[index]=$("root:reference"+num2str(index))
	endfor
	variable/G V_FitTol = V_FitTol_input
	//the fit itself
	for (row_index=v_rowstart;row_index<v_rowend;row_index+=1)
		for (column_index=0;column_index<v_columns_max;column_index+=1)
 			make/D/O/N=(v_num_of_waves) W_coef 	
 			killwaves W_coef
 			duplicate W_coef_input W_coef
 			make/D/O/N=(v_num_of_waves) W_sigma=NaN
			make/D/N=(v_layers_max)/O tempwave=0
			make/D/N=(v_layers_max)/O tempwave2=0
			MatrixOP/O/Z tempwave=beam(w_sample_duplicate,row_index,column_index)
			setscale/P x,dimoffset(w_sample_duplicate,2),dimdelta(w_sample_duplicate,2),tempwave,tempwave2
			FuncFit/Q/H=s_hold_LCM   f_fitFunc_lc_MT_AllAtOnce_singleshift W_coef tempwave[v_start,v_end] /M=w_mask /NWOK // /H=s_hold_LCM 
			//save coefficients
			for (v_num_index=0;v_num_index<(v_num_of_waves);v_num_index+=1)
				w_export[row_index][column_index][v_num_index]=W_coef[v_num_index]
				w_stddev[row_index][column_index][v_num_index]=W_sigma[v_num_index]
			endfor
		endfor
	endfor
end

function f_ThreadGroupRelease_all()
	variable dummy
	dummy= ThreadGroupRelease(-2)
end

//Contraints text to thread safe waves
function f_constr_txt_to_matrix_MT(string importstring, variable number_of_Ks) //Generates input for multi thread funcfit
	importstring=replacestring("K",importstring,"")
	importstring=replacestring("\"",importstring,"")
	importstring=replacestring(",",importstring,";")
	variable v_items=itemsinList(importstring)
	make/O/N=(v_items) temporary_wave
	variable index
	for (index=0;index<v_items;index+=1)
		temporary_wave[index]=str2num(StringFromList(index,importstring))
	endfor
	make/O/N=(v_items,number_of_Ks) M_FitConstraint=0
	make/O/N=(v_items) W_FitConstraint
	string tempstring	
	string s_replace
	variable v_row=0
	variable v_read
	for (index=0;index<v_items;index+=1)
		tempstring=StringFromList(index,importstring)
		if (StringMatch(tempstring,"*>*")==1)
			s_replace=num2str(temporary_wave[index])+" >"
			tempstring=replacestring(s_replace,tempstring,"")
			v_read=str2num(tempstring)
			M_FitConstraint[v_row][temporary_wave[index]]=-1
			if (v_read==0)
				W_FitConstraint[v_row]=v_read
			else
				W_FitConstraint[v_row]=(-1)*v_read
			endif
		else
			s_replace=num2str(temporary_wave[index])+" <"
			tempstring=replacestring(s_replace,tempstring,"")
			v_read=str2num(tempstring)
			M_FitConstraint[v_row][temporary_wave[index]]=1
			W_FitConstraint[v_row]=v_read
		endif
		v_row+=1
	endfor
	killwaves temporary_wave	
end

function f_change_LCM_start_values(variable v_input) //change all non-offset initial guesses
	string s_control=ControlNameList("LCMpanel", ";", "V_start_*")
	variable index
	for (index=0;index<(ItemsInList(s_control,";")-1);index+=1)	
		SetVariable $(StringFromList(index,s_control,";")),win=LCMpanel ,value=_NUM:v_input
	endfor
end

function f_change_LCM_min_values(variable v_input) //change all non-offset initial guesses
	string s_control=ControlNameList("LCMpanel", ";", "V_Min*")
	variable index
	for (index=0;index<(ItemsInList(s_control,";")-1);index+=1)	
		SetVariable $(StringFromList(index,s_control,";")),win=LCMpanel ,value=_NUM:v_input
	endfor
end

function f_change_LCM_min_values_all(variable v_input) //change all non-offset initial guesses
	string s_control=ControlNameList("LCMpanel", ";", "V_Min*")
	variable index
	for (index=0;index<(ItemsInList(s_control,";"));index+=1)	
		SetVariable $(StringFromList(index,s_control,";")),win=LCMpanel ,value=_NUM:v_input
	endfor
end

function f_change_LCM_max_values(variable v_input) //change all non-offset initial guesses
	string s_control=ControlNameList("LCMpanel", ";", "V_Max*")
	variable index
	for (index=0;index<(ItemsInList(s_control,";"));index+=1)	
		SetVariable $(StringFromList(index,s_control,";")),win=LCMpanel ,value=_NUM:v_input
	endfor
end


//the fit function
//normal LCM
Threadsafe function f_fitFunc_lc_MT(w,x) : FitFunc
    Wave w
    Variable x
    wave /WAVE list=root:listwave
    int index, imax
    variable v_return=0
 	 imax = dimsize(list,0)
    for (index=0;index<imax;index++)
   		wave ref_part=list[index]
    	v_return+=w[index]*ref_part(x)
    endfor
    v_return+=w[(dimsize(list,0))]
    return v_return
end


Threadsafe function f_fitFunc_lc_MT_AllAtOnce(wave pw, wave yw, wave xw) : FitFunc
    wave /WAVE list=root:listwave
    int index, imax
    imax = dimsize(list,0)
    yw=0
    for (index=0;index<imax;index++)
        wave w = list[index]
        yw += pw[index] * w(xw[p])
    endfor
    yw += pw[index]
end

//multiple shift LCM
Threadsafe function f_fitFunc_lc_MT_AllAtOnce_multishift(wave pw, wave yw, wave xw) : FitFunc
    wave /WAVE list=root:listwave
    int index, imax
    imax = dimsize(list,0)
    yw=0
    for (index=0;index<imax;index++)
        wave w = list[index]
        yw += pw[index] * w(xw[p]+pw[(imax+1+index)])
    endfor
    yw += pw[index]
end


Threadsafe function f_fitFunc_lc_MT_multishift(w,x) : FitFunc
    Wave w
    Variable x
    wave /WAVE list=root:listwave
    int index, imax
    variable v_return=0
 	 imax = dimsize(list,0)
    for (index=0;index<imax;index++)
   		wave ref_part=list[index]
    	v_return+=w[index]*ref_part(x+w[imax+1+index])
    endfor
    v_return+=w[(dimsize(list,0))]
    return v_return
end

//single shift LCM
Threadsafe function f_fitFunc_lc_MT_singleshift(w,x) : FitFunc
    Wave w
    Variable x
    wave /WAVE list=root:listwave
    int index, imax
    variable v_return=0
 	 imax = dimsize(list,0)
    for (index=0;index<imax;index++)
   		wave ref_part=list[index]
    	v_return+=w[index]*ref_part(x+w[imax+1])
    endfor
    v_return+=w[(dimsize(list,0))]
    return v_return
end

Threadsafe function f_fitFunc_lc_MT_AllAtOnce_singleshift(wave pw, wave yw, wave xw) : FitFunc
    wave /WAVE list=root:listwave
    int index, imax
    imax = dimsize(list,0)
    yw=0
    for (index=0;index<imax;index++)
        wave w = list[index]
        yw += pw[index] * w(xw[p]+pw[(imax+1)])
    endfor
    yw += pw[index]
end

function f_sigma_to_sigma_rel()
	string s_sigma=removeEnding(getbrowserSelection(0),"sigma")+"sigma"
 	if (waveexists($s_sigma)==0)
 		print "Error: No sigma_wave selected"
 		return 0
 	endif
 	wave w_sigma=$s_sigma
 	string s_coef=removeEnding(s_sigma,"sigma")+"coef"
  	if (waveexists($s_coef)==0)
 		print "Error: No coef_wave selected"
 		return 0
 	endif	
 	wave w_coef=$s_coef
 	duplicate/FREE w_sigma w_sigma_rel
 	w_sigma_rel=w_sigma/w_coef
 	string s_sigma_rel=s_sigma+"_rel"
 	duplicate/O w_sigma_rel $s_sigma_rel
end