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