1subroutine lorentzian(y,npts,a)
2
3! Input:  y(npts); assume x(i)=i, i=1,npts
4! Output: a(5)
5!         a(1) = baseline
6!         a(2) = amplitude
7!         a(3) = x0
8!         a(4) = width
9!         a(5) = chisqr
10
11  real y(npts)
12  real a(5)
13  real deltaa(4)
14
15  a=0.
16  df=12000.0/8192.0                               !df = 1.465 Hz
17  width=0.
18  ipk=0
19  ymax=-1.e30
20  do i=1,npts
21     if(y(i).gt.ymax) then
22        ymax=y(i)
23        ipk=i
24     endif
25!     write(50,3001) i,i*df,y(i)
26!3001 format(i6,2f12.3)
27  enddo
28!  base=(sum(y(ipk-149:ipk-50)) + sum(y(ipk+51:ipk+150)))/200.0
29  base=(sum(y(1:20)) + sum(y(npts-19:npts)))/40.0
30  stest=ymax - 0.5*(ymax-base)
31  ssum=y(ipk)
32  do i=1,50
33     if(ipk+i.gt.npts) exit
34     if(y(ipk+i).lt.stest) exit
35     ssum=ssum + y(ipk+i)
36  enddo
37  do i=1,50
38     if(ipk-i.lt.1) exit
39     if(y(ipk-i).lt.stest) exit
40     ssum=ssum + y(ipk-i)
41  enddo
42  ww=ssum/y(ipk)
43  width=2
44  t=ww*ww - 5.67
45  if(t.gt.0.0) width=sqrt(t)
46  a(1)=base
47  a(2)=ymax-base
48  a(3)=ipk
49  a(4)=width
50
51! Now find Lorentzian parameters
52
53  deltaa(1)=0.1
54  deltaa(2)=0.1
55  deltaa(3)=1.0
56  deltaa(4)=1.0
57  nterms=4
58
59!  Start the iteration
60  chisqr=0.
61  chisqr0=1.e6
62  do iter=1,5
63     do j=1,nterms
64        chisq1=fchisq0(y,npts,a)
65        fn=0.
66        delta=deltaa(j)
6710      a(j)=a(j)+delta
68        chisq2=fchisq0(y,npts,a)
69        if(chisq2.eq.chisq1) go to 10
70        if(chisq2.gt.chisq1) then
71           delta=-delta                      !Reverse direction
72           a(j)=a(j)+delta
73           tmp=chisq1
74           chisq1=chisq2
75           chisq2=tmp
76        endif
7720      fn=fn+1.0
78        a(j)=a(j)+delta
79        chisq3=fchisq0(y,npts,a)
80        if(chisq3.lt.chisq2) then
81           chisq1=chisq2
82           chisq2=chisq3
83           go to 20
84        endif
85
86! Find minimum of parabola defined by last three points
87        delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5)
88        a(j)=a(j)-delta
89        deltaa(j)=deltaa(j)*fn/3.
90!          write(*,4000) iter,j,a,chisq2
91!4000      format(i1,i2,4f10.4,f11.3)
92     enddo
93     chisqr=fchisq0(y,npts,a)
94!       write(*,4000) 0,0,a,chisqr
95     if(chisqr/chisqr0.gt.0.99) exit
96     chisqr0=chisqr
97  enddo
98  a(5)=chisqr
99
100  return
101end subroutine lorentzian
102
103