c
c     COPYRIGHT (c) 2008
c     
c     Applied Research Laboratories, The University of Texas at Austin
c     (ARL:UT) All rights reserved.
c 
c     This material may be reproduced by or for the U.S. Government
c     only. This notice must appear in all copies of this file and its
c     derivatives. Further distribution outside the U.S. Government is
c     prohibited without the explicit written authorization from the
c     Director of the Environmental Sciences Laboratory, ARL:UT.
c 
c     THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
c     EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE
c     IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
c     PURPOSE, OR NONINFRINGEMENT.
c 
c     IN NO EVENT WILL ARL:UT BE LIABLE FOR ANY DAMAGES WHATSOEVER
c     ARISING OUT OF THE USE OR INABILITY TO USE THIS SOFTWARE, EVEN
c     IF ARL:UT HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
c
subroutine data_cfar_sdn_mex( &
	indata, &	! 1 real*8 (ntime)
	ntime, &		! 2 integer*4
	nt_low, &		! 3 integer*4
	nt_high, &		! 4 integer*4
	nt_guard1, &	! 5 integer*4
	nt_guard2, &	! 6 integer*4
	pct_cens, &		! 7 real*8
	mean_fac, &		! 8 real*8
	SDN, &			! 9 real*8 (ndB,nnu)
	ndB, &			! 10 integer*4
	nnu, &			! 11 integer*4
	dB0, &			! 12 real*8
	ddB, &			! 13 real*8
	nu0, &			! 14 real*8
	dnu, &			! 15 real*8
	outdata_dB)		! 16 real*8 (ntime)

! LAB cfar a matrix

implicit none

integer (kind=4) :: ntime,dim2,nt_low,nt_high,nt_guard1,nt_guard2,ndB,nnu
real (kind=8), dimension(ntime) :: indata, outdata_dB
real (kind=8) :: pct_cens,mean_fac,dB0,ddB,nu0,dnu,drmean,drstd
real (kind=8), dimension(ndB,nnu) :: SDN

integer (kind=4) :: jtime,ntmin,ntmax,npt_avg,lii,i,lmu,jnu,jdB
integer (kind=4), dimension(ntime) :: ii
real (kind=8), dimension(ntime) :: data_sort,windata,mudata
real (kind=8) :: data_avg,mu,sigma,nu,dval



!do j2=1,dim2

	! un-dB
!	do jtime=1,ntime
!		indata(jtime)=10.**(indata_dB(jtime)/10.)
!	enddo

	do jtime=1,ntime

		! get various window indices
	    ntmin=max(1,jtime-nt_low)
	    ntmax=ntmin+nt_low+nt_high
	    if (ntmax>ntime) then
	    	ntmax=ntime 
	    	ntmin=max(1,ntmax-nt_low-nt_high)
	    endif
		!ii=[ntmin:jtime-nt_guard1-1,jtime+nt_guard2+1:ntmax];
		lii=0
		do i=ntmin,jtime-nt_guard1-1
			lii=lii+1
			ii(lii)=i
		enddo
		do i=jtime+nt_guard2+1,ntmax
			lii=lii+1
			ii(lii)=i
		enddo
	    npt_avg=int((1 - .01*pct_cens)*lii)
		!select data in the window
		do i=1,lii
			windata(i)=indata(ii(i))
			data_sort(i)=windata(i) ! make a copy for hpsort to clobber
		enddo
!        [data_sort]=sort(indata);
!    	call hpsort(lii,data_sort(1:lii))
    	call hpsort(lii,data_sort)
!        data_avg=mean(data_sort(1:npt_avg));
		data_avg=drmean(data_sort(1:npt_avg),npt_avg)
!        indx_use=find(indata<mean_fac*data_avg);
!        mu=mean(indata(indx_use)); 
		lmu=0
		do i=1,lii
			if (windata(i)<mean_fac*data_avg) then
				lmu=lmu+1
				mudata(lmu)=windata(i)
			endif
		enddo
		mu=drmean(mudata(1:lmu),lmu)
!        if mu~=1, sigma=std(indata(indx_use)); else, sigma=1; end
		sigma=drstd(mudata(1:lmu),lmu) ! ignore above test for now;  watch for editted input = 0 dB => 1
!        nu=max(1,min(20,(mu/sigma)^2));
		nu=max(1.,min(20.,real((mu/sigma)**2)))
!        KD{kdim}=jtime;
!        dval=10*log10(indata0(KD{1},KD{2},KD{3},KD{4},KD{5},KD{6},KD{7})/mu);
		dval=10*log10(indata(jtime)/mu)
!        jnu=round((nu-nu0)/dnu)+1;
		jnu=nint((nu-nu0)/dnu)+1
!        jdB=min(ndB,max(1,round((dval-dB0)/ddB)+1));
		jdB=min(ndB,max(1,nint((dval-dB0)/ddB)+1))
!        %%disp([jdB,jnu]);
!        SP.data(KD{1},KD{2},KD{3},KD{4},KD{5},KD{6},KD{7})=SDN.data(jdB,jnu);
		outdata_dB(jtime)=SDN(jdB,jnu)


!write(4,*) 'jtime   data_avg   mu    sigma   nu    dval   jnu   jdB'
!write(4,*) jtime, data_avg, mu, sigma, nu, dval, jnu, jdB
!write(4,*) 'indata'
!write(4,'(20(f6.2,1x))') indata
!write(4,*) 'windata'
!write(4,'(20(f6.2,1x))') windata(1:lii)
!write(4,*) 'data_sort'
!write(4,'(20(f6.2,1x))') data_sort(1:lii)

	enddo
!enddo


return
end



      subroutine hpsort(n,ra)
!
      implicit none
      integer*4 n
      integer*4 i,ir,j,l
      real*8 ra(n)
!: Sorts an array ra(1:n) into ascending order using the Heapsort
!: algorithm. n is input; ra is replaced on output by its sorted 
!: rearrangement.
!
      real*8 rra
      if(n .lt. 2) return
!: The index l will be decremented from its initial value down to 1
!: during the "hiring" (heap creation) phase.  Once it reaches 1, the
!: index ir will be decremented from its initial value down to 1
!: during the "retirement and promotion" (heap selection) phase.
      l=n/2 + 1
      ir=n
10    continue
         if(l .gt. 1) then
            l=l-1
            rra=ra(l)
         else
            rra=ra(ir)
            ra(ir)=ra(1)
            ir=ir-1
            if(ir .eq. 1) then
               ra(1)=rra
               return
            endif
         endif
         i=l
         j=l+l
20       if(j .le. ir) then
            if(j .lt. ir) then
               if(ra(j) .lt. ra(j+1)) j=j+1
            endif
            if(rra .lt. ra(j)) then
               ra(i)=ra(j)
               i=j
               j=j+j
            else
               j=ir+1
            endif
            goto 20
         endif
         ra(i)=rra
      goto 10

      end



function drmean(dat,n)
! double real mean

integer (kind=4) :: n,i
real (kind=8), dimension(n) :: dat
real (kind=8) :: drmean


drmean=0
do i=1,n
	drmean=drmean+dat(i)
enddo
drmean=drmean/float(n)

return
end function drmean



function drstd(dat,n)
!double real unbias standard deviation

integer (kind=4) :: n,i
real (kind=8), dimension(n) :: dat
real (kind=8) :: drstd,sum,mean


mean=drmean(dat,n)
sum=0
do i=1,n
	sum=sum+(dat(i)-mean)**2
enddo
drstd=sqrt(sum/(n-1))

return
end




      STDCALL SUBROUTINE MEXFUNCTION(NLHS, PLHS, NRHS, PRHS)

      use ABF90

      integer plhs(*), prhs(*)
!     argument pointers
      integer indata_dB,ntime,dim2,nt_low,nt_high,nt_guard1,nt_guard2, &
         pct_cens,mean_fac,SDN,ndB,nnu,dB0,ddB,nu0,dnu,outdata_dB

      integer nlhs, nrhs
      integer m, n
!     Check for proper number of arguments. 
      if(nrhs .ne. 15) then
         call mexErrMsgTxt('15 inputs required.')
      elseif(nlhs .ne. 1) then
         call mexErrMsgTxt('One output required.')
      endif

!     Check to ensure the array is numeric (not strings).
      if(mxIsNumeric(prhs(1)) .eq. 0) then
         call mexErrMsgTxt('Input 1 must be a double vector.')
      endif
      if(mxIsNumeric(prhs(2)) .eq. 0) then
         call mexErrMsgTxt('Input 2 must be a int32.')
      endif
      if(mxIsNumeric(prhs(3)) .eq. 0) then
         call mexErrMsgTxt('Input 3 must be a int32.')
      endif
      if(mxIsNumeric(prhs(4)) .eq. 0) then
         call mexErrMsgTxt('Input 4 must be a int32.')
      endif
      if(mxIsNumeric(prhs(5)) .eq. 0) then
         call mexErrMsgTxt('Input 5 must be a int32.')
      endif
      if(mxIsNumeric(prhs(6)) .eq. 0) then
         call mexErrMsgTxt('Input 6 must be a int32.')
      endif
      if(mxIsNumeric(prhs(7)) .eq. 0) then
         call mexErrMsgTxt('Input 7 must be a double real.')
      endif
      if(mxIsNumeric(prhs(8)) .eq. 0) then
         call mexErrMsgTxt('Input 8 must be a double real.')
      endif
      if(mxIsNumeric(prhs(9)) .eq. 0) then
         call mexErrMsgTxt('Input 9 must be a double array.')
      endif
     if(mxIsNumeric(prhs(10)) .eq. 0) then
         call mexErrMsgTxt('Input 10 must be a int32.')
      endif
     if(mxIsNumeric(prhs(11)) .eq. 0) then
         call mexErrMsgTxt('Input 11 must be a int32.')
      endif
     if(mxIsNumeric(prhs(12)) .eq. 0) then
         call mexErrMsgTxt('Input 12 must be a double real.')
      endif
     if(mxIsNumeric(prhs(13)) .eq. 0) then
         call mexErrMsgTxt('Input 13 must be a double real.')
      endif
     if(mxIsNumeric(prhs(14)) .eq. 0) then
         call mexErrMsgTxt('Input 14 must be a double real.')
      endif
     if(mxIsNumeric(prhs(15)) .eq. 0) then
         call mexErrMsgTxt('Input 15 must be a double real.')
      endif

!     assign pointers for input arguments
	indata_dB = mxGetPr(prhs(1))
	ntime = mxGetPr(prhs(2))
	nt_low = mxGetPr(prhs(3))
	nt_high = mxGetPr(prhs(4))
	nt_guard1 = mxGetPr(prhs(5))
	nt_guard2 = mxGetPr(prhs(6))
	pct_cens = mxGetPr(prhs(7))
	mean_fac = mxGetPr(prhs(8))
	SDN = mxGetPr(prhs(9))
	ndB = mxGetPr(prhs(10))
	nnu = mxGetPr(prhs(11))
	dB0 = mxGetPr(prhs(12))
	ddB = mxGetPr(prhs(13))
	nu0 = mxGetPr(prhs(14))
	dnu = mxGetPr(prhs(15))

!     Get the size of the input array.
      m = mxGetM(prhs(1))
      n = mxGetN(prhs(1))
!     Create matrix for the return argument.
      plhs(1) = mxCreateFull(m, n, 0)
      outdata_db = mxGetPr(plhs(1))

!     Call the computational subroutine.
	call data_cfar_sdn_mex( &
	%val(indata_dB), &	! 1 real*8 (ntime)
	%val(ntime), &		! 2 integer*4
	%val(nt_low), &		! 4 integer*4
	%val(nt_high), &		! 5 integer*4
	%val(nt_guard1), &	! 6 integer*4
	%val(nt_guard2), &	! 7 integer*4
	%val(pct_cens), &		! 8 real*8
	%val(mean_fac), &		! 9 real*8
	%val(SDN), &			! 10 real*8 (ndB,nnu)
	%val(ndB), &			! 11 integer*4
	%val(nnu), &			! 12 integer*4
	%val(dB0), &			! 13 real*8
	%val(ddB), &			! 14 real*8
	%val(nu0), &			! 15 real*8
	%val(dnu), &			! 16 real*8
	%val(outdata_dB) )		! 17 real*8 (ntime)

      return
      end
