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