diff --git a/lib/determ.f90 b/lib/determ.f90 new file mode 100644 index 000000000..6ae471bd3 --- /dev/null +++ b/lib/determ.f90 @@ -0,0 +1,32 @@ +real*8 function determ(array,norder) + implicit real*8 (a-h,o-z) + real*8 array(10,10) + + determ=1. + do k=1,norder + if (array(k,k).ne.0) go to 41 + do j=k,norder + if(array(k,j).ne.0) go to 31 + enddo + determ=0. + go to 60 + +31 do i=k,norder + s8=array(i,j) + array(i,j)=array(i,k) + array(i,k)=s8 + enddo + determ=-1.*determ +41 determ=determ*array(k,k) + if(k.lt.norder) then + k1=k+1 + do i=k1,norder + do j=k1,norder + array(i,j)=array(i,j)-array(i,k)*array(k,j)/array(k,k) + enddo + enddo + end if + enddo + +60 return +end function determ diff --git a/lib/flat3.f90 b/lib/flat3.f90 new file mode 100644 index 000000000..8255fca5c --- /dev/null +++ b/lib/flat3.f90 @@ -0,0 +1,74 @@ +subroutine flat3(savg0,iz,nterms,ynoise,savg) + + implicit real*8 (a-h,o-z) + parameter (NSMAX=6827) + real*4 savg0(iz) + real*4 savg(iz) + real*4 ynoise,y4 + + real*8 x(NSMAX) + real*8 y(NSMAX) + real*8 y0(NSMAX) + real*8 yfit(NSMAX) + real*8 a(10) + integer ii(NSMAX) + + npts0=999999 + df=12000.0/16384.0 + + do i=1,iz + y0(i)=db(savg0(i)) + enddo + ia=200.0/df + ib=4500.0/df + j=0 + do i=ia,ib + j=j+1 + x(j)=i*df + y(j)=y0(i) + ii(j)=i + enddo + + npts=j + mode=0 + a=0.0 + + do iter=1,99 + call polfit(x,y,y,npts,nterms,mode,a,chisqr) +! print*,iter,npts,a(1:nterms) + + rewind 21 + do i=1,ib + f=i*df + yfit(i)=0.0 + do n=1,nterms + yfit(i)=yfit(i) + a(n) * f**(n-1) + enddo +! write(21,1010) f,y0(i),yfit(i),y0(i)-yfit(i) +!1010 format(4f12.3) + enddo + k=0 + do j=1,npts + y1=y(j)-yfit(ii(j)) + if(y1.lt.ynoise) then + k=k+1 + x(k)=x(j) + y(k)=y(j) + ii(k)=ii(j) + endif + enddo + npts=k + if(npts.eq.npts0 .or. npts.lt.(ib-ia)/10) exit + npts0=npts + enddo + +! do j=1,npts +! write(22,1010) x(j),y(j),yfit(ii(j)),y(j)-yfit(ii(j)) +! enddo + + do i=1,ib + y4=y0(i)-yfit(i) + savg(i)=10.0**(0.1*y4) + enddo + +end subroutine flat3 diff --git a/lib/polfit.f90 b/lib/polfit.f90 new file mode 100644 index 000000000..ea7ffc9f6 --- /dev/null +++ b/lib/polfit.f90 @@ -0,0 +1,72 @@ +subroutine polfit(x,y,sigmay,npts,nterms,mode,a,chisqr) + implicit real*8 (a-h,o-z) + real*8 x(npts), y(npts), sigmay(npts), a(nterms) + real*8 sumx(19), sumy(10), array(10,10) + +! Accumulate weighted sums + nmax = 2*nterms-1 + sumx=0. + sumy=0. + chisq=0. + do i=1,npts + xi=x(i) + yi=y(i) + if(mode.lt.0) then + weight=1./abs(yi) + else if(mode.eq.0) then + weight=1 + else + weight=1./sigmay(i)**2 + end if + xterm=weight + do n=1,nmax + sumx(n)=sumx(n)+xterm + xterm=xterm*xi + enddo + yterm=weight*yi + do n=1,nterms + sumy(n)=sumy(n)+yterm + yterm=yterm*xi + enddo + chisq=chisq+weight*yi**2 + enddo + +! Construct matrices and calculate coefficients + do j=1,nterms + do k=1,nterms + n=j+k-1 + array(j,k)=sumx(n) + enddo + enddo + + delta=determ(array,nterms) + if(delta.eq.0) then + chisqr=0. + a=0. + else + do l=1,nterms + do j=1,nterms + do k=1,nterms + n=j+k-1 + array(j,k)=sumx(n) + enddo + array(j,l)=sumy(j) + enddo + a(l)=determ(array,nterms)/delta + enddo + +! Calculate chi square + + do j=1,nterms + chisq=chisq-2*a(j)*sumy(j) + do k=1,nterms + n=j+k-1 + chisq=chisq+a(j)*a(k)*sumx(n) + enddo + enddo + free=npts-nterms + chisqr=chisq/free + end if + + return +end subroutine polfit diff --git a/lib/symspec.f90 b/lib/symspec.f90 index 2facb4903..a08e3c1a7 100644 --- a/lib/symspec.f90 +++ b/lib/symspec.f90 @@ -40,7 +40,8 @@ subroutine symspec(k,ntrperiod,nsps,ingain,slope,pxdb,s,df3,ihsym,npts8) go to 900 !Wait for enough samples to start endif - if(nfft3.ne.nfft3z .or. slope.ne.slope0) then !New nfft3, compute window + if(nfft3.ne.nfft3z .or. (slope.ne.slope0 .and. abs(slope+0.1).gt.0.05)) then +! Compute new window and adjust scale factor pi=4.0*atan(1.0) do i=1,nfft3 w3(i)=2.0*(sin(i*pi/nfft3))**2 !Window for nfft3 spectrum @@ -103,10 +104,14 @@ subroutine symspec(k,ntrperiod,nsps,ingain,slope,pxdb,s,df3,ihsym,npts8) s(i)=gain*sx enddo -! s=0.05*s/ref s=scale*s savg=scale*ssum/ihsym + if(abs(slope+0.1).lt.0.01) then + call flat3(s,iz,3,1.0,s) + call flat3(savg,iz,3,1.0,savg) + endif + 900 npts8=k/8 return diff --git a/wsjtx.pro b/wsjtx.pro index 69ff9ece2..418c2e5d7 100644 --- a/wsjtx.pro +++ b/wsjtx.pro @@ -7,7 +7,7 @@ QT += network multimedia greaterThan(QT_MAJOR_VERSION, 4): QT += widgets CONFIG += thread -CONFIG += console +#CONFIG += console TARGET = wsjtx #DESTDIR = ../qt4_install