-- AutoDeconv.lua 	automatic deconvolution of a kinetic experiment
-- this script has been written for a specific spectrum, yet it can be easily adapted to similar experiments

-- DEFINITIONS

HOME = HOME or "/Users/GB/"		 -- define here your home directory
						-- never forget the initial slash !
filename = "dec_table_auto.txt"	--here we save the table of integrals

NumSpectra = 32	-- number of points for the kinetic study
NumRegions = 5	-- regions to deconvolute

Fi = {}		-- central frequencies at time zero
Ff = {}		-- central frequencies at the end of the study
par = {}	-- parameter for the deconvolution
-- they must be obtained from a preliminary deconvolution performed on the first spectrum
HalfSpan = {}	-- region to deconvolute, in ppm units, based upon experience and observation
-- the regions we'll submit to deconvolution will be 2 * HalfSpan wide
-- the value depends on the shape of the signal (e.g.: a multiplet requires a larger region than a singlet)

local i = 1	-- progressive index, simplifies the editing of this script
-- for example, you can reorder the definitions below and they will still work

-- proton 4
Fi[i] = 3.7104
Ff[i] = 3.7114
HalfSpan[i] = 0.020
par[i] = [[
Parameters for 1 peaks
  frequency (Hz)     intensity      width (Hz)   Lorentzian %

      1855.4433       150.0000         1.9975       100.0000
]]
i = i + 1

-- proton 5
Fi[i] = 3.9368
Ff[i] = 3.9185
HalfSpan[i] = 0.026
par[i] = [[
Parameters for 2 peaks
  frequency (Hz)     intensity      width (Hz)   Lorentzian %

      1963.1406        0.5000         4.7583       100.0000
      1955.7102        0.5000         4.8810       100.0000
]]
i = i + 1

-- proton 6
Fi[i] = 5.6611
Ff[i] = 5.6734
HalfSpan[i] = 0.020
par[i] = [[
Parameters for 1 peaks
  frequency (Hz)     intensity      width (Hz)   Lorentzian %

      2833.0753         2.0000         2.5857       100.0000
]]
i = i + 1

-- proton 7
Fi[i] = 5.0177
Ff[i] = 4.9926
HalfSpan[i] = 0.024
par[i] = [[
Parameters for 2 peaks
  frequency (Hz)     intensity      width (Hz)   Lorentzian %

      2512.4897         1.0000         4.0698       100.0000
      2505.4187         1.0000         3.7240       100.0000
]]
i = i + 1

-- proton 8
Fi[i] = 3.75
Ff[i] = 3.76
HalfSpan[i] = 0.025
par[i] = [[
Parameters for 1 peaks
  frequency (Hz)     intensity      width (Hz)   Lorentzian %

      1874.5991         8.000         3.6089       100.0000
]]
i = i + 1

-- end of definitions ----------------------------------------

-- ************->    HERE WE GO:    <-***************

io.output(HOME..filename)		-- create/open the file were the results will be stored
io.write("point")	-- header
for i = 1,NumRegions do
	io.write("\tppm\t\tarea")
end
io.write("\n")

spectral = getf("x")				-- read some experimental parameters 
conversion = 1.0 / spectral.MHz		-- useful to convert from Hz to ppm
spectral = getf("y")
step = spectral.width / spectral.size
Y = spectral.start + step * (spectral.size - 0.01)	-- position of the first row on its dummy ppm scale

for P = 1,NumSpectra do			-- NumSpectra
	io.write( string.format("\n\t%02d\t", P ) )	-- report the experiment no.
	mark('h', Y )		-- choose a row
	extract()			-- extract the coresponding 1D spectrum
	delint()			-- we need to normalize the intensities
	region( 2.4, 2.6 )	-- region containing protons no. 1 and 9
	press 'i'			-- first integral, automatically set to 1
	intreg( 1, 300 )	-- we set it to 300 to have manageable numbers (>1 and <300)
	for R = 1,NumRegions do
		local C = Fi[R] + (Ff[R] - Fi[R]) * (P-1)/(NumSpectra-1)
		region( C + HalfSpan[R], C - HalfSpan[R] )
		cliptext( par[R] )
		local onv = dec()
		local NumPeaks = 0
		if onv then
			dec( onv, "paste" )		-- input our starting parameters
			dec( onv, "sall" )		-- fit all parameters
			dec( onv, "SAME" )		-- Lorentzian line-shapes
			dec( onv, "ok  " )		-- Fit
			dec( onv, "ok  " )		-- refine
			dec( onv, "copy" )		-- output refined parameters
			NumPeaks = dec( onv )
			dec( onv, "close" )
		end
		local str = cliptext()
		local area = 0
		local frequency = 0
		local pos = string.find( str, "\n\n" )
		local mistake = not pos
		if pos then
			pos = pos + 1
			for i = 1,NumPeaks do
				pos = string.find( str, "\n", pos )
				mistake = not pos
				if mistake then print 'no endline' break end
				pos = pos + 1
				local F
				local A
				_, _, F, A = string.find( str, "([%w.]+)%s+([%w.]+)", pos )
				mistake = not F or not A
				if mistake then print 'no values' break end
				frequency = frequency + F
				area = area + A
			end
			if NumPeaks > 0 then
				frequency = conversion * frequency / NumPeaks
			end
		end
		if mistake then		-- information for debugging
			print("Error Processing Region no.", R, "of Spectrum no.", P)
			print("deconvolution returned:")
			print( str )
		end
		if R == 1 or P > 2 then
			par[R] = str
		end
		io.write( string.format("%0.3f\t%0.3f\t", frequency, area ) )
	end
	closex()
	Y = Y - step		-- position of the next spectrum
	print("Point", P, "completed")
end

io.close()