1subroutine polfit(y,npts,a)
2
3! Input:  y(npts)                !Expect npts=4
4! Output: a(1) = baseline
5!         a(2) = amplitude
6!         a(3) = theta (deg)
7
8  real y(npts)
9  real a(3)
10  real deltaa(3)
11  integer ipk(1)
12  save
13
14! Set starting values:
15  a(1)=minval(y)
16  a(2)=maxval(y)-a(1)
17  ipk=maxloc(y)
18  a(3)=(ipk(1)-1)*45.0
19
20  deltaa(1:2)=0.1*a(2)
21  deltaa(3)=10.0
22  nterms=3
23
24!  Start the iteration
25  chisqr=0.
26  chisqr0=1.e6
27  iters=10
28
29  do iter=1,iters
30     do j=1,nterms
31        chisq1=fchisq_pol(y,npts,a)
32        fn=0.
33        delta=deltaa(j)
3410      a(j)=a(j)+delta
35        chisq2=fchisq_pol(y,npts,a)
36        if(chisq2.eq.chisq1) go to 10
37        if(chisq2.gt.chisq1) then
38           delta=-delta                      !Reverse direction
39           a(j)=a(j)+delta
40           tmp=chisq1
41           chisq1=chisq2
42           chisq2=tmp
43        endif
4420      fn=fn+1.0
45        a(j)=a(j)+delta
46        chisq3=fchisq_pol(y,npts,a)
47        if(chisq3.lt.chisq2) then
48           chisq1=chisq2
49           chisq2=chisq3
50           go to 20
51        endif
52
53! Find minimum of parabola defined by last three points
54        delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5)
55        a(j)=a(j)-delta
56        deltaa(j)=deltaa(j)*fn/3.
57!          write(*,4000) iter,j,a,deltaa,chisq2
58!4000      format(2i2,2(2x,3f8.2),f12.5)
59     enddo  ! j=1,nterms
60     chisqr=fchisq_pol(y,npts,a)
61!     write(*,4000) 0,0,a,chisqr
62     if(chisqr.lt.1.0) exit
63     if(deltaa(1).lt.0.01*(a(2)-a(1)) .and. deltaa(2).lt.0.01*(a(2)-a(1))   &
64          .and. deltaa(3).lt.1.0) exit
65     if(chisqr/chisqr0.gt.0.99) exit
66     chisqr0=chisqr
67  enddo  ! iter
68  a(3)=mod(a(3)+360.0,180.0)
69
70  return
71end subroutine polfit
72
73real function fchisq_pol(y,npts,a)
74
75  real y(npts),a(3)
76  data rad/57.2957795/
77
78  chisq = 0.
79  do i=1,npts
80     theta=(i-1)*45.0
81     yfit=a(1) + a(2)*cos((theta-a(3))/rad)**2
82     chisq=chisq + (y(i) - yfit)**2
83  enddo
84  fchisq_pol=chisq
85
86  return
87end function fchisq_pol
88