C+
C NAME:
C	NODAT
C PURPOSE:
C	NOrmalize DATa: subtract a zodiacal dust cloud model from the raw data
C CATEGORY:
C	Data processing
C CALLING SEQUENCE:
	program NoDat
C INPUTS:
C	Unnormalized Helios data file (.ZLD file).
C OUTPUTS:
C	Output file has the same name as input file, except for the file type
C	which is .DAT and is created in the same directory as the .ZLD file.
C CALLS:
C	AskYN, DustAsymmetry, LsqFit, bOpenFile, iHOSRead, iHOSWrite, HOSRead
C	HOSUpdate, Say, ArrI4AddConstant, BadR4, iHOSInfo, iFreeLun
C	iPutFileSpec, iGetFileSpec, iGetSymbol, iSetSymbol, iDeleteSymbol
C INCLUDE:
	include		'filparts.h'
	include		'openfile.h'
	include		'hos_e9.h'
C RESTRICTIONS:
C	The maximum number of points fitted is given by the parameter value nT
C	(1000 at present). The polarization brightness is NOT fitted.
C PROCEDURE:
C	LsqFit at end of this file
C >	The zodiacal dust cloud brightness for each line of sight 
C	(photometer/sector) combination is modeled by the function
C	I = C*r^P (I = intensity [S10], r= s/c-Sun distance [AU]). The
C	constants C and P least squares fitted (after taking
C	logarithms, the fit becomes linear).
C >	After the fit is removed, the average intensity over the residual time
C	series for each photometer/sector/filter combination is subtracted 
C	(i.e. the 'background' level is redefined)
C >	The color index is reset from 101,102,103 in the .ZLD input file to
C	1,2,3 in the .DAT output file
C MODIFICATION HISTORY:
C	1990, Valerie Hung. Modification of earlier version with the same name.
C	JUN-1992, Paul Hick (UCSD); Added option to use correction for 
C			asymmetry of zodiacal dust cloud
C-
	parameter	(nT=1000)
	parameter	(nS=32)
	integer		N(nT)
	real		T(nT)
	real		R(nT)
	real		L(nT)
	integer		P(nT)
	integer		C(nT)
	integer		F(nT)
	real		Z(nS,nT)

	integer		Lst(3,11)
	integer		nSAll(3)	/nS,nS,2/

	integer		kPP(3)		/HOS__P_1,HOS__P_2,HOS__P_3/
	integer		kCC(3)		/HOS__C_1,HOS__C_2,HOS__C_3/
	integer		kFF(5)		/HOS__F_1,HOS__F_2,HOS__F_3,HOS__F_4,HOS__F_5/

	character	cFile1*80
	character	cFile2*80
	character	cSay*5		/'NoDat'/
	character	cSymbol*17	/'LIB__SIGNAL_QUIET'/
	character	cQuiet*5

	logical		bOpenFile
	logical		bSym

	Bad = BadR4()

	iR = iHOSInfo(0,cFile1,iRecl,iSc,iYr,TMin,TMax,kPCF,Lst)
	if (iR .eq. 0) call Say(cSay,'S','StopPlay','Bye bye')
	if (iand(kPCF,HOS__NONORM) .eq. 0) call Say(cSay,'E','File','already normalized')

	I = iPutFileSpec(FIL__TYPE,FIL__TYPE,'.dat')
	I = iGetFileSpec(0,FIL__TYPE,cFile2)

	bSym = bOpenFile(OPN__HOS+OPN__REOPEN+OPN__STOP+OPN__READONLY,iU1,cFile1,iRecl)
	bSym = bOpenFile(OPN__HOS+OPN__REOPEN+OPN__STOP+OPN__NEW     ,iU2,cFile2,iRecl)

	iR = 0
	do while (iHOSRead(kPCF,iU1,iRecl,iR,T,P,C,F,R,L,nS,Z) .eq. 0)
	    I = iHOSWrite(kPCF,iU2,iRecl,iR,T,P,C,F,R,L,nS,Z)
	end do

	iU1 = iFreeLun(iU1)
	iU2 = iFreeLun(iU2)

	call AskYN('Correct for asymmetry in dust cloud$yes',bSym)

	iQ = iGetSymbol(cSymbol,cQuiet)
	I  = iSetSymbol(cSymbol,'I',2)		! Suppress informational messages

	! Loop over 33 photom/color/filter combinations
	! 16 deg photometer: 5 filters (skip pB), three colors: 15 combinations
	! 31 deg photometer: 5 filters (skip pB), three colors: 15 combinations
	! 90 deg photometer: 2 filter, 3 colors: 3 combinations

	do kPF=1,11				! Photometers and filters
	    kP = 1+(kPF-1)/5
	    kF = 4
	    if (kP .ne. 3) kF = 1+mod(kPF-1,5)

	    do kC=1,3				! Colors
		if (Lst(kC,kPF) .ne. 0) then	! # Records
		    write (*,'(3(3X,A,I2))') 'Photometer:',kP,'Color:',kC,'Filter:',kF

		    kPCF = kPP(kP)+kCC(kC)+kFF(kF)
		    call HOSRead(kPCF,cFile2,nT,iT,N,T,P,C,F,R,L,nS,Z)
						! Set color to 1,2,3
		    call ArrI4AddConstant(iT,C,-100,C)
		    do I=1,iT
			N(I) = IBSET(N(I),31)	! Flag records as modified
		    end do
						! Correct for assymetry in dust cloud
		    if (bSym) call DustAsymmetry(0,iSc,iT,P,L,nS,Z)

		    do iS=1,nSAll(kP)		! Loop over sectors
			if (kF .ne. 5 .and. .not. (kP .eq. 3 .and. iS .eq. 2)) call LsqFit(iS,iT,R,nS,Z)
						! Fit and subtract power law (unless pB)
			nP = 0
			Zt = 0.
			do I=1,iT
			    if (Z(iS,I) .ne. Bad) then
				nP = nP+1
				Zt = Zt+Z(iS,I)
			    end if
			end do
			if (nP .ne. 0) then	! Subtract average
			    Zt = Zt/nP
			    do I=1,iT
				if (Z(iS,I) .ne. Bad) Z(iS,I) = Z(iS,I)-Zt
			    end do
			end if
		    end do

		    call HOSUpdate(kPCF-HOS__NONORM,cFile2,nT,iT,N,T,P,C,F,R,L,nS,Z)
		end if
	    end do
	end do

	I = iDeleteSymbol(cSymbol,2)
	if (iQ .ne. 0) iQ = iSetSymbol(cSymbol,cQuiet,2)

	call Say(cSay,'S','StopPlay','Normalized file is '//cFile2)

	end

	subroutine LsqFit(iS,iT,R,nS,Z)		! Fit Z ~ A*R^G
	    integer	iS
	    integer	iT
	    real	R(iT)
	    integer	nS
	    real	Z(nS,iT)		! or log(Z) ~ log(A)+G*log(R)

	Bad = BadR4()

	sR  = 0.
	sZ  = 0.
	sRZ = 0.
	sRR = 0.
	do I=1,iT
	    if (Z(iS,I) .ne. Bad) then
		aR  = alog(R(I))
		aZ  = alog(Z(iS,I))
		sR  = sR +aR
		sZ  = sZ +aZ
		sRZ = sRZ+aR*aZ
		sRR = sRR+aR*aR
	    end if
	end do

	A =   iT*sRR-sR*sR
	G = ( iT*sRZ-sZ*sR )/A
	A = (sRR*sZ -sR*sRZ)/A
	A = exp(A)

	do I=1,iT			! Subtract power law fit
	    if (Z(iS,I) .ne. Bad) Z(iS,I) = Z(iS,I)-A*R(I)**G
	end do

	return
	end