```
; please separate the individual routines----------------------------------

function one_gs,x,y,a,da
;
; function to fit one Gaussian with zero point (baseline)
; usage: yfit=two_gs(x,y,a,da)
; input: x,y,a
; output: yfit,da
; da(4) are diff.corr.'s to parameters a(4) as follows:
;   a(0) - central strength
;   a(1) - central position
;   a(2) - width
;   a(3) - zero point

; a(0) and a(3) - central strengths
; a(1) and a(4) - central positions
; a(2) and a(5) - widths
; a(6) is zero point
;

n = n_elements(x)
c = x & e1=x
c = ((x-a(1))/a(2))^2
e1(*)=0.
d = where(c le 30.)
e1(d) = exp(-c(d))
res = y - a(0)*e1 - a(3)
print,'Sum res^2 = ',total(res^2)
;

t = fltarr(n,4)         ; derivatives
t(*,0) = e1
t(*,1) = 2*e1*a(0)*(x-a(1))/a(2)^2
t(*,2) = t(*,1)*(x-a(1))/a(2)
t(*,3) = 1.
;
tt = transpose(t)#t     ; solution
tr = transpose(t)#res
svd,tt,w,u,v
wp = fltarr(4,4)
mw = max(w)
for i=0,3 do if w(i) ge mw*1.e-6 then wp(i,i)=1./w(i)
da = v#wp#(transpose(u)#tr)
a1 = a + da
;
c = x & e1=x            ; new Gaussian
e1(*)=0.
c = ((x-a1(1))/a1(2))^2
d = where(c le 30.)
e1(d) = exp(-c(d))
res = y - a1(0)*e1 - a1(3)
print,'Sum new res^2 = ',total(res^2)
return,a1(0)*e1+a1(3)
end
; end of one_gs.pro

; -------------------------------------------------------------------------
function two_gs,x,y,a,da
;
; function to fit two Gaussians
; usage: yfit=two_gs(x,y,a,da)
; input: x,y,a
; output: yfit,da
; da are diff.corr.'s to parameters a
; a(0) and a(3) - central strengths
; a(1) and a(4) - central positions
; a(2) and a(5) - widths
; a(6) is zero point
;
n = n_elements(x)
c = x & e1=x & e2=x
c = ((x-a(1))/a(2))^2
e1(*)=0. & e2(*)=0.
d = where(c le 30.)
e1(d) = exp(-c(d))
c = ((x-a(4))/a(5))^2
d = where(c le 30.)
e2(d) = exp(-c(d))
res = y - a(0)*e1 - a(3)*e2-a(6)
print,'Sum res^2 = ',total(res^2)
;
t = fltarr(n,7)
t(*,0) = e1
t(*,3) = e2
t(*,1) = 2*e1*a(0)*(x-a(1))/a(2)^2
t(*,4) = 2*e2*a(3)*(x-a(4))/a(5)^2
t(*,2) = t(*,1)*(x-a(1))/a(2)
t(*,5) = t(*,4)*(x-a(4))/a(5)
t(*,6) = 1.
;
tt = transpose(t)#t
tr = transpose(t)#res
svd,tt,w,u,v
wp = fltarr(7,7)
mw = max(w)
for i=0,6 do if w(i) ge mw*1.e-6 then wp(i,i)=1./w(i)
da = v#wp#(transpose(u)#tr)
a1 = a + da
;
c = x & e1=x & e2=x
e1(*)=0. & e2(*)=0.
c = ((x-a1(1))/a1(2))^2
d = where(c le 30.)
e1(d) = exp(-c(d))
c = ((x-a1(4))/a1(5))^2
d = where(c le 30.)
e2(d) = exp(-c(d))
res = y - a1(0)*e1 - a1(3)*e2-a1(6)
print,'Sum new res^2 = ',total(res^2)
return,a1(0)*e1+a1(3)*e2+a1(6)
end
; end of two_gs.pro

;------------------------------------------------------------------------------
function two_gs1,x,y,a,da
;
; function to fit two Gaussians
; usage: yfit=two_gs1(x,y,a,da)
; input: x,y,a
; output: yfit,da
; da are diff.corr.'s to parameters a
; a(0) and a(3) - central strengths
; a(1) and a(4) - central positions
; a(2) and a(5) - widths                 [a(5) kept constant]
; a(6) is zero point
;
; NOTE: this version keeps the width of the 2nd Gaussian fixed
; but inputs and outputs are the same as for "two_gs"
;
n = n_elements(x)
c = x & e1=x & e2=x
c = ((x-a(1))/a(2))^2
e1(*)=0. & e2(*)=0.
d = where(c le 30.)
e1(d) = exp(-c(d))
c = ((x-a(4))/a(5))^2
d = where(c le 30.)
e2(d) = exp(-c(d))
res = y - a(0)*e1 - a(3)*e2-a(6)
print,'Sum res^2 = ',total(res^2)
;
t = fltarr(n,6)
t(*,0) = e1
t(*,3) = e2
t(*,1) = 2*e1*a(0)*(x-a(1))/a(2)^2
t(*,4) = 2*e2*a(3)*(x-a(4))/a(5)^2
t(*,2) = t(*,1)*(x-a(1))/a(2)   ; to solve for width of 1st Gaussian
t(*,5) = 1.
;
tt = transpose(t)#t
tr = transpose(t)#res
svd,tt,w,u,v
wp = fltarr(6,6)
mw = max(w)
for i=0,5 do if w(i) ge mw*1.e-6 then wp(i,i)=1./w(i)
dda = v#wp#(transpose(u)#tr)
da = [dda(0:4),0.,dda(5)]   ; width of 2nd Gaussian unchanged
a1 = a + da
;
c = x & e1=x & e2=x
e1(*)=0. & e2(*)=0.
c = ((x-a1(1))/a1(2))^2
d = where(c le 30.)
e1(d) = exp(-c(d))
c = ((x-a1(4))/a1(5))^2
d = where(c le 30.)
e2(d) = exp(-c(d))
res = y - a1(0)*e1 - a1(3)*e2-a1(6)
print,'Sum new res^2 = ',total(res^2)
return,a1(0)*e1+a1(3)*e2+a1(6)
end

; end of two_gs1.pro
; -----------------------------------------------------------------------
function two_gs2,x,y,a,da
;
; function to fit two Gaussians
; usage: yfit=two_gs2(x,y,a,da)
; input: x,y,a (values of x,y for function and 7 parameters a)
; output: yfit,da (fit and corrections to a)
; a(0) and a(3) - central strengths
; a(1) and a(4) - central positions
; a(2) and a(5) - widths       [both kept constant]
; a(6) is zero point
;
; NOTE: this version keeps widths of both Gaussians fixed
;    as full solutions normally unstable

n = n_elements(x)
c = x & e1=x & e2=x
e1(*)=0. & e2(*)=0.
c = ((x-a(1))/a(2))^2
d = where(c le 30.)
e1(d) = exp(-c(d))
c = ((x-a(4))/a(5))^2
d = where(c le 30.)
e2(d) = exp(-c(d))
res = y - a(0)*e1 - a(3)*e2-a(6)
print,'Sum res^2 = ',total(res^2)
;
t = fltarr(n,5)  ; widths not determined
t(*,0) = e1
t(*,2) = e2
t(*,1) = 2*e1*a(0)*(x-a(1))/a(2)^2
t(*,3) = 2*e2*a(3)*(x-a(4))/a(5)^2
t(*,4) = 1.
;
tt = transpose(t)#t   ; solution
tr = transpose(t)#res
svd,tt,w,u,v
wp = fltarr(5,5)
mw = max(w)
for i=0,4 do if w(i) ge mw*1.e-6 then wp(i,i)=1./w(i)
dda = v#wp#(transpose(u)#tr)          ; corrections computed
da = [dda(0:1),0.,dda(2:3),0.,dda(4)] ; no changes to widths
a1 = a + da
;
c = x & e1=x & e2=x  ; initialize
e1(*)=0. & e2(*)=0.
c = ((x-a1(1))/a1(2))^2
d = where(c le 30.)
e1(d) = exp(-c(d))   ; first Gaussian
c = ((x-a1(4))/a1(5))^2
d = where(c le 30.)
e2(d) = exp(-c(d))   ; second Gaussian
res = y - a1(0)*e1 - a1(3)*e2-a1(6)  ; residuals
print,'Sum new res^2 = ',total(res^2)
return,a1(0)*e1+a1(3)*e2+a1(6)
end
; end of two_gs2.pro and whole package
;---------------------------------------------------------------------
```