rm(list=ls())
parameter.filename = "Rinput.csv";
new.parameter.filename = "Routput.csv";

params = scan(file=parameter.filename,skip=1,what=list(1.0,1.0,1.0,1.0,""),sep=",",strip.white=T)

a = as.numeric(params[1])
b = as.numeric(params[2])
t0 = as.numeric(params[3])
V = as.numeric(params[4])
quadratic.peaks.filename = as.character(params[5])


peakdata = read.csv(file=quadratic.peaks.filename, header = TRUE, sep = ",", quote="\"", dec=".", fill = TRUE)


# target mass values are mis (perhaps from reference spectrum)
# actual mass values are pis
mis = peakdata[,1]
pis = peakdata[,2]

# try to make sure only valid numeric values are input
mis = mis[is.numeric(mis) & !is.na(mis)]
pis = pis[is.numeric(pis) & !is.na(pis)]


# function to minimize
leastsqssum3pscale = function(pars,masses,tofs,V){
  a = pars[1]; b = pars[2]; t0 = pars[3];
  leastsqsum3p = sum( ((masses - V*(a*((tofs - t0)^2)+b))/masses)^2)
  leastsqsum3p
}

# calculate the tis in the paper, here called tofpis
# masses too near to zero can cause errors if they correspond
# to negative times of flight -- in that case need to use
# parts of the spectrum away from 0
tofpis = sqrt((pis/V-b)/a)+t0
pars3 = c(a,b,t0); vals = pars3

# find minimizing values of a, b, t0
vals = optim(pars3,leastsqssum3pscale,masses=mis,tofs=tofpis,V=V)$par

anew = vals[1]; bnew = vals[2]; t0new = vals[3]

# new m/z values are given by the next line
#Mz = V*anew*(tofpis - t0new)^2 + V*bnew

# write out the data
newdf = data.frame(Avalue = anew, Bvalue = bnew, t0value = t0new, Uvalue = V);

write.table(file=new.parameter.filename,newdf,sep=",",
            col.names = TRUE,row.names=F,quote=F)

q("no")