C	DCJG.FOR - to generate time ordered MAP data
C	and analyze it without using data when reference beam is in the
C	galactic plane
C
	PARAMETER (JRES=8,NP=12*4**JRES,NBIG=4*NP)
	PARAMETER (DTR = 57.2957795130824)
C	only data with reference beam at sin|b| > SINBLIM is used
	REAL*4 SINBLIM
	LOGICAL*1 OK(NP)
	LOGICAL OLD
	COMMON /MASK/ OK
	INTEGER*2 COVER(NP)
	REAL*8 BB
	REAL*8 VXM(NP),GUESS(NP),RHS(NP)
	REAL*8 ESTIMATE(NP),RESID(NP),PP(NP)
	REAL*4 DGMX,EMAX0,EMAX
	REAL*4 GALAXY(NBIG)
C	GUESS and VXM are not needed while creating COVER and RHS
C	so they can be stored in the space used for the simulated map
	EQUIVALENCE (GALAXY,VXM),(GALAXY(2*NP+1),RESID)
C	VXM and GUESS are not needed simultaneously
	EQUIVALENCE (VXM,GUESS)
	REAL*8 T68,T,TS,GMT2T68
	REAL*8 RHO,RHO_OLD,PQ,BETA,RHOVERPQ
	REAL*4 WALL,WALLS
	COMPLEX*16 QSIMG,QG
	PARAMETER (NGMAX=200)
	INTEGER JPIX(NGMAX)
	REAL*4 B(4,NGMAX),A(NGMAX)
	REAL*4 C(4,NGMAX)
	REAL*4 GS(3,10),G(3),GD(3)
	LOGICAL RESTART
	CHARACTER*1 ANS
	CHARACTER*80 FN
	WRITE (*,898) 'Is this a restart? (Y/N):'
	READ (*,899) NC,ANS
	RESTART = ANS.EQ.'Y' .OR. ANS.EQ.'y'
C	read the simulated map
	WRITE (*,898) 'Enter res=9 simulated map file:'
898	FORMAT(1X,A,$)
	READ (*,899) NC,FN
899	FORMAT(Q,A)
	IF (.NOT.RESTART) THEN
	  OPEN (UNIT=1,FILE=FN(1:NC),STATUS='OLD',READONLY,
	1	FORM='UNFORMATTED')
	  READ (1) (GALAXY(I),I=1,NBIG)
	  CLOSE (UNIT=1)
	ENDIF
C
	WRITE (*,898) 'Enter noise/sample & IX,IY,IZ:'
	READ (*,FMT='(E10.3,3I6)') SGMA,KX,KY,KZ
	WRITE (*,898) 'Enter list of beam positions:'
	READ (*,899) NC,FN
	OPEN (UNIT=2,FILE=FN(1:NC),STATUS='OLD',READONLY)
C	main + beam first
	READ (2,900) (B(I,1),I=1,3)
	B(4,1) = 0
	CALL NORM(B)
C	copy + beam with y coordinate reversed to make - beam
	DO I=1,4
	  B(I,2) = B(I,1)
	ENDDO
	B(2,2) = -B(2,2)
	A(1) = 1
	A(2) = -1
900	FORMAT(3F10.5,2E11.3)
C	read list of near sidelobes
	NG=2
10	READ (2,900,END=20) (B(I,NG+1),I=1,3),A(NG+1)
	  NG=NG+1
	  CALL NORM(B(1,NG))
	  B(4,NG) = 0
C	copy the plus side for the minus side
	  DO I=1,4
	    B(I,NG+1) = B(I,NG)
	  ENDDO
C	reverse the sign
	  A(NG+1) = -A(NG)
	  NG = NG+1
C	reverse the y coordinate
	  B(2,NG) = -B(2,NG)
	  GO TO 10
20	CLOSE (UNIT=2)
	WRITE (*,*) 'Beams:'
	DO J=1,NG
	  WRITE (*,FMT='(3F10.6,1PE11.3)') (B(I,J),I=1,3),A(J)
	ENDDO
	WRITE (*,898) 'Enter starting time (YYDDDHHMMSSFFF):'
	READ (*,899) NC,FN
	FN = FN(1:NC)//'00000000000000'
	TS = GMT2T68(FN(1:14))
	WRITE (*,898) 'Enter Nday, hours:'
	READ (*,*) NDAY,NHOUR
	WRITE (*,*) 'Simulating ',NDAY,' days.  Every ',NHOUR,'th hour.'
	IF (RESTART) GO TO 50
	WRITE (*,898) 'Enter sin|b| min:'
	READ (*,*) SINBLIM
	NS = 0
30	WRITE (*,898) 'Enter (l,b) of source to exclude:'
	READ (*,*,END=40) GLON,GLAT
	NS = NS+1
	GS(3,NS) = SIN(GLAT/DTR)
	GS(1,NS) = COS(GLAT/DTR)*COS(GLON/DTR)
	GS(2,NS) = COS(GLAT/DTR)*SIN(GLON/DTR)
	GO TO 30
40	WRITE (*,*) '# of sources to exclude:',NS
C	the formula below gives negative values for negative inputs
C	negative values mean NO exclusion
	DVSLIM = SINBLIM*ABS(SINBLIM)
C	set up reference beam exclusion map
	DO I=1,NP
	  CALL NESTED_CENPIX(I-1,JRES,G)
C	set OK to true outside plane and false in the galactic plane
	  OK(I) = ABS(G(3)).GT.SINBLIM
C	set OK to false near the listed sources
	  DO J=1,NS
	    IF (DVS(G,GS(1,J)).LT.DVSLIM) OK(I) = .FALSE.
	  ENDDO
	ENDDO
C
C	do the simulation
C
50	WRITE (*,*) 'Tstart=',TS
	WALLS = SECNDS(0.0)
C	on restart just read in the RHS and COVER
C	ignore the "no path" compiler warning
	IF (RESTART) THEN
	  READ (2) COVER
	  READ (3) RHS
	  READ (4) OK
	  CLOSE (UNIT=2)
	  CLOSE (UNIT=3)
	  CLOSE (UNIT=4)
	  READ (11,ERR=60,END=60) (GALAXY(I),I=1,NP)
	  DO I=1,NP
	    ESTIMATE(I) = GALAXY(I)
	    IF (.NOT.OK(I)) ESTIMATE(I)=0
	  ENDDO
	  CLOSE (UNIT=11)
	  GO TO 100
C	checkpoint not found, use RHS/diag(A) as guess
60	  RESTART = .FALSE.
	  GO TO 100
	ENDIF
C	MAKE_RHS generates "real" [simulated] data and makes the RHS and COVER
C	vectors
70	CALL MAKE_RHS(TS,NDAY,NHOUR,2,NG,B,A,GALAXY,OK,SGMA,COVER,RHS)
	DT = SECNDS(WALLS)
	WRITE (*,*) 'COVER & RHS done at',DT,' seconds.'
C	write out the results
	WRITE (2) COVER
	CLOSE (UNIT=2)
	WRITE (3) RHS
	CLOSE (UNIT=3)
	WRITE (4) OK
	CLOSE (UNIT=4)
C
C	initialize map making
C
100	BB = 0
	DO I=1,NP
	  BB = BB+RHS(I)**2
C	extend exclusion zone to the un-covered parts of the map
	  IF (COVER(I).LE.0.5) OK(I) = .FALSE.
	ENDDO
	BB = SQRT(BB)
	WRITE (*,*) 'L2 norm of whole RHS = ',BB
C	first guess for solution of M*X = RHS is ESTIMATE = RHS/diag(M)
	IF (.NOT.RESTART) CALL V_DIV_DIAG(COVER,RHS,ESTIMATE)
C	get initial error for convergence testing
	CALL V_NORM(RHS,EMAX0,OK)
	CALL V_NORM(ESTIMATE,DGMX,OK)
	DT = SECNDS(WALLS)
	WRITE (*,*) DT,EMAX0,DGMX
C	find RESID = RHS - M*ESTIMATE
C	AXV multiplies matrix M times vector ESTIMATE giving vector RESID
	CALL AXV(TS,NDAY,NHOUR,2,B,A,OK,RESID,ESTIMATE)
	DO I=1,NP
	  RESID(I) = RHS(I)-RESID(I)
	ENDDO
	CALL V_NORM(RESID,EMAX,OK)
	CALL V_DIV_DIAG(COVER,RESID,GUESS)
	CALL V_NORM(GUESS,DGMX,OK)
	DT = SECNDS(WALLS)
	WRITE (*,*) DT,EMAX,DGMX
C	on first iteration of conjugate gradient, use beta=0
190	BETA = 0
	OLD = .FALSE.
C	top of iteration loop for conjugate gradient
200	RHO = 0
	DO I=1,NP
C	note we are only using the off-plane "OK" part of the vectors
	  IF(OK(I)) RHO = RHO+RESID(I)*GUESS(I)
	ENDDO
	IF (OLD) BETA = RHO/RHO_OLD
C
C	with exclusion zone the matrix structure is
C
C	/       \ /   \   /   \
C	| A   0 | | x |   | b |
C	|       | |   | = |   |
C	| C   D | | y |   | e |
C	\       / \   /   \   /
C
C	where A is symmetric and D is diagonal.  x is the part of the map 
C	that is outside the galaxy exclusion zone, while y is the galaxy
C	exclusion zone.
C
C	We first solve Ax = b by a conjugate gradient technique:
C
C	x = b/diag(A)
C	r = b-Ax
C	g = r/diag(A)
C
C	while ||r|| too big do
C	  rho = r*g
C	  if first beta = 0 else beta = rho/rho'
C	  p = g + beta*p
C	  q = Ap
C	  x = x+[rho/(p*q)]*p
C	  r = r-[rho/(p*q)]*p
C	  rho' = rho
C	endwhile
C
C	We then fill in the galaxy part of the map using
C	y = (E-Cx)/D
C
250	WRITE (*,*) 'beta,rho=',BETA,RHO
	DO I=1,NP
	  PP(I) = GUESS(I) + BETA*PP(I)
	ENDDO
	CALL AXV(TS,NDAY,NHOUR,2,B,A,OK,VXM,PP)
	PQ = 0
	DO I=1,NP
	  IF (OK(I)) PQ = PQ+VXM(I)*PP(I)
	ENDDO
	CALL V_NORM(PP,DGMX,OK)
	RHOVERPQ = RHO/PQ
	DGMX =  RHOVERPQ*DGMX
	DO I=1,NP
	  ESTIMATE(I) = ESTIMATE(I) + RHOVERPQ*PP(I)
	  RESID(I) = RESID(I) - RHOVERPQ*VXM(I)
	ENDDO
	WRITE (*,*) 'rho/PQ, PQ:',RHOVERPQ,PQ
	CALL V_NORM(RESID,EMAX,OK)
	DT = SECNDS(WALLS)
	WRITE (*,*) DT,EMAX,DGMX
C
C	checkpoint write of ESTIMATE, but use R*4 copy
C	can use GALAXY which is equivalenced to VXM and GUESS - none are needed
	DO I=1,NP
	  IF (OK(I)) THEN
	    GALAXY(I) = ESTIMATE(I)
	  ELSE
	    GALAXY(I) = -1.E30
	  ENDIF
	ENDDO
	OPEN (UNIT=11,STATUS='UNKNOWN',FORM='UNFORMATTED')
	WRITE (11) (GALAXY(I),I=1,NP)
	CLOSE (UNIT=11)
C
C	test for convergence
C
	IF (EMAX.LT.3.162E-8*EMAX0) GO TO 300
290	RHO_OLD = RHO
	OLD = .TRUE.
	CALL V_DIV_DIAG(COVER,RESID,GUESS)
	GO TO 200
C
C	ending.  recompute the residuals
300	WRITE (*,*) 'checking residuals again'
	DO I=1,NP
	  IF (.NOT.OK(I)) ESTIMATE(I) = 0
	ENDDO
	CALL AXV(TS,NDAY,NHOUR,2,B,A,OK,VXM,ESTIMATE)
	DO I=1,NP
	  RESID(I) = RHS(I)-VXM(I)
	ENDDO
	CALL V_NORM(RESID,EMAX,OK)
	IF (EMAX.GT.5.0E-8*EMAX0) GO TO 290
	WRITE (*,*) 'filling in the Galaxy'
	DO I=1,NP
	  IF (COVER(I).GT.0.5 .AND. .NOT.OK(I)) THEN
	    ESTIMATE(I) = RESID(I)/COVER(I)
	  ENDIF
	  IF (COVER(I).LT.0.5) ESTIMATE(I) = -1.0E30
	  GALAXY(I) = ESTIMATE(I)
	ENDDO
	WRITE (11) (GALAXY(I),I=1,NP)
	STOP
	END

	SUBROUTINE MAKE_RHS(TS,NDAY,NHOUR,NSLV,NSM,B,A,GALAXY,OK,
	1	SIGMA,COVER,RHS)
C
C	computes the time-ordered data for a MAP simulation and generates
C	RHS and COVER
C	Inputs:	TS, the time of first data point, in seconds since 5/24/68
C	NDAY: the number of days to simulate minus 1: loop runs from 0 to NDAY
C	NHOUR: spacing in hours between hours simulated.  should evenly
C		divide 24.  NHOUR=1 means do every hours, 2 means do every
C		other hour, 3 means do 8 hours per day.  This variable is used
C		to reduce computation time during experimental runs
C	NSLV:	not used
C	NSM:	total number of beams, main+sidelobes
C	B:	positions of the beams
C	A:	amplitudes of beams.  A(1) and A(2) must be +1 and -1
C	GALAXY:	map with 3145728 HEALpix pixels of the input sky
C	OK:	galaxy exclusion mask
C	SIGMA:	RMS noise per sample
C
C	Outputs: COVER, the number of times each pixel is observed
C	RHS:	the "right hand side" of the equation M*X = RHS to be solved
C		for the map
C
C	external routines used:
C	QSIMG(T):	returns a quaternion describing the spacecraft attitude
C		at time T in seconds after 24 May 1968.  Note the FORTRAN
C		is told that quaterions are COMPLEX*16 quantities
C	DIPOLE(T,GD):	returns the dipole vector in mK in (l,b) coordinates
C		at time T
C	BODY2CEL(QG,NSM,B,C): rotates the NSM positions in array B by the
C		rotation specified by the quaterion QG giving a result in
C		array C
C	NESTED_PIXNO(C,IRES,JP): converts position given by unit vector C
C		into a HEALpixe pixel number at resolution IRES, or
C		for n = 12*4^IRES pixels.  JP is output between 0 and n-1
C	
	PARAMETER (JRES=8,NP=12*4**JRES,NBIG=4*NP)
	PARAMETER (DTR = 57.2957795130824)
	LOGICAL*1 OK(NP)
	REAL*8 TS,T
	INTEGER*2 COVER(NP)
	REAL*8 RHS(NP)
	REAL*4 GALAXY(NBIG)
	PARAMETER (NGMAX=200)
	INTEGER JPIX(NGMAX)
	REAL*4 B(4,NGMAX),A(NGMAX)
	REAL*4 C(4,NGMAX)
	REAL*4 GD(3)
	COMPLEX*16 QG,QSIMG
	DO I=1,NP
	  COVER(I) = 0
	ENDDO
	DO I=1,NP
	  RHS(I) = 0
	ENDDO
C
	DO ID=0,NDAY
C	simulating every NHOUR'th hour
	 DO IH=0,23,NHOUR
	  T = TS+86400*ID+3600*IH+1800
	  CALL DIPOLE(T,GD)
	  DO IS=0,3599
	    DO JF=0,4
	      T = TS+86400*ID+3600*IH+IS+JF/5.D0
	      QG = QSIMG(T)
C	plus beam
	      CALL BODY2CEL(QG,NSM,B,C)
	      CALL NESTED_PIXNO(C,9,JP9)
	      JPLUS = JP9/4
C	minus beam
	      CALL NESTED_PIXNO(C(1,2),9,JM9)
	      VD = DOT(C,GD)-DOT(C(1,2),GD)
	      JMINUS = JM9/4
C	generate the signals, sum the COUNTs and RHS's
C	main beam signal
	      VM = GALAXY(JP9+1)-GALAXY(JM9+1)
	      VSL = 0.
	      DO I=3,NSM
	        CALL NESTED_PIXNO(C(1,I),9,JPIX(I))
	        VSL = VSL+A(I)*GALAXY(JPIX(I)+1)
	      ENDDO
	      VM = VM+VSL+VD+SIGMA*GAUSRN(IR,JR)
	      IF (OK(JMINUS+1)) THEN
	        COVER(JPLUS+1) = COVER(JPLUS+1)+1
	        RHS(JPLUS+1) = RHS(JPLUS+1)+VM
	      ENDIF
	      IF (OK(JPLUS+1)) THEN
	        COVER(JMINUS+1) = COVER(JMINUS+1)+1
	        RHS(JMINUS+1) = RHS(JMINUS+1)-VM
	      ENDIF
	    ENDDO	! JF
	  ENDDO		! IS
	 ENDDO		! IH
	ENDDO		! ID
	RETURN
	END

	SUBROUTINE AXV(TS,NDAY,NHOUR,NSLV,B,A,OK,VXM,V)
C
C	multiplies the vector V by the "matrix" giving vector VXM
C	V is used as a map, the simulation is flown over it, and VXM
C	is what would have been RHS in the routine MAKE_RHS
C
C	Inputs:	TS, the time of first data point, in seconds since 5/24/68
C	NDAY: the number of days to simulate minus 1: loop runs from 0 to NDAY
C	NHOUR: spacing in hours between hours simulated.  should evenly
C		divide 24.  NHOUR=1 means do every hours, 2 means do every
C		other hour, 3 means do 8 hours per day.  This variable is used
C		to reduce computation time during experimental runs
C	NSLV:	not used
C	B:	positions of the beams
C	A:	amplitudes of beams.  A(1) and A(2) must be +1 and -1
C	OK:	galaxy exclusion mask
C
C	external routines used:
C	QSIMG(T):	returns a quaternion describing the spacecraft attitude
C		at time T in seconds after 24 May 1968.  Note the FORTRAN
C		is told that quaterions are COMPLEX*16 quantities
C	BODY2CEL(QG,NSM,B,C): rotates the NSM positions in array B by the
C		rotation specified by the quaterion QG giving a result in
C		array C
C	NESTED_PIXNO(C,IRES,JP): converts position given by unit vector C
C		into a HEALpixe pixel number at resolution IRES, or
C		for n = 12*4^IRES pixels.  JP is output between 0 and n-1
C	
	PARAMETER (JRES=8,NP=12*4**JRES,NBIG=4*NP)
	PARAMETER (DTR = 57.2957795130824)
	LOGICAL*1 OK(NP)
	REAL*8 TS,T
	REAL*8 VXM(NP),V(NP)
	PARAMETER (NGMAX=200)
	INTEGER JPIX(NGMAX)
	REAL*4 B(4,NGMAX),A(NGMAX)
	REAL*4 C(4,NGMAX)
	COMPLEX*16 QG,QSIMG
	DO I=1,NP
	  VXM(I) = 0
	ENDDO
C
	DO ID=0,NDAY
C	simulating every NHOUR'th hour
	 DO IH=0,23,NHOUR
	  DO IS=0,3599
	    DO JF=0,4
	      T = TS+86400*ID+3600*IH+IS+JF/5.D0
	      QG = QSIMG(T)
C	plus beam
	      CALL BODY2CEL(QG,NSLV,B,C)
	      CALL NESTED_PIXNO(C,JRES,JPLUS)
C	minus beam
	      CALL NESTED_PIXNO(C(1,2),JRES,JMINUS)
	      IF (OK(JMINUS+1)) THEN
	        VXM(JPLUS+1) = VXM(JPLUS+1) + (V(JPLUS+1)-V(JMINUS+1))
	      ENDIF
	      IF (OK(JPLUS+1)) THEN
	        VXM(JMINUS+1) = VXM(JMINUS+1)-(V(JPLUS+1)-V(JMINUS+1))
	      ENDIF
	    ENDDO	! JF
	  ENDDO		! IS
	 ENDDO		! IH
	ENDDO		! ID
	RETURN
	END

	SUBROUTINE V_DIV_DIAG(COVER,RHS,ESTIMATE)
	PARAMETER (JRES=8,NP=12*4**JRES)
	INTEGER*2 COVER(NP)
	REAL*8 RHS(NP),ESTIMATE(NP)
	REAL*8 SUM,WSUM
	SUM = 0
	WSUM = 0
	DO I=1,NP             
	  IF (COVER(I).GT.0.5) THEN
	    ESTIMATE(I) = RHS(I)/COVER(I)
	    WSUM = WSUM + COVER(I)
	    SUM = SUM+COVER(I)*ESTIMATE(I)
	  ELSE
	    ESTIMATE(I) = 0
	  ENDIF
	ENDDO
	SUM = SUM/WSUM
	DO I=1,NP             
	  IF (COVER(I).GT.0.5) THEN
	    ESTIMATE(I) = ESTIMATE(I) - SUM
	  ENDIF
	ENDDO
	RETURN
	END

	SUBROUTINE V_NORM(V,E,OK)
	PARAMETER (JRES=8,NP=12*4**JRES)
	REAL*8 V(NP),SUM
	REAL*4 E
	LOGICAL*1 OK(NP)
	SUM=0
	DO I=1,NP
	  IF (OK(I)) SUM=SUM+V(I)**2
	ENDDO
	E = SQRT(SUM)
	RETURN
	END

	SUBROUTINE DIPOLE(T68,G)
C	returns the dipole vector in galactic coordinates for time T68
C	based on motion on the COBE DMR 4 year results:
C	dipole = 3.353(24) mK, l=264.25(33), b=48.22(13)
	REAL*8 T68	! input, seconds since 5/24/1968
	REAL*4 G(3)	! output, units are milliKelvin
	REAL*8 T0
	PARAMETER (T0 = 1020384.D3)	! T68 of year 2000 autumnal equinox
C	PARAMETER (BETAT0 = 0.2725)	! (v/c)*T_\circ of CMB
	PARAMETER (BETAT0 = 0.0)	! TESTing TESTing TESTing
	REAL*8 EARTHLONG		! heliocentric long of Earth in radians
	REAL*4 E(3),C(3),GV(3)
C	at autumnal equinox, Sun is at 180^o longitude or 12h RA so
C	Earth is at 0^o long and moving in the +y direction
C	90 days later Earth is at 90^o and moving in -x direction
	EARTHLONG = 1.991062161890534D-7*(T68-T0)
	E(3) = 0
	E(2) = BETAT0*COS(EARTHLONG)
	E(1) = -BETAT0*SIN(EARTHLONG)
C	tan(long) for E is -cos(el)/sin(el) so long of E is EARTHLONG+pi/2
	CALL ETOC(E,C)
	CALL CTOG(C,GV)
C	(lambda,beta) of COBE dipole: (171.8,-11.2)
	G(3) = GV(3)+2.5004
	G(2) = GV(2)-2.2228
	G(1) = GV(1)-0.2238
C	max mag of G occurs 83 days after autumnal equinox = 12/15/2000
	RETURN
	END
