;+
;
;  Name:
;       DET_POINT
;
;
; PURPOSE:
;  	Compute direction cosines for BATSE LAD's for Solar Viewing
;	FIND POINTING CHANGES IN THE SPACECRAFT X-AXIS
;
;
; CATEGORY: cgro/batse, attitude analysis
;
;
; CALLING SEQUENCE:
;	det_point, ut, x_ra, x_dec, pnt_ndx, cos_det,z_ra,z_dec,$
;	  celestial_fov=batsefov, cos_spec=cos_spec, $
;	  spacecraftfov=batse, x2fov=x2fov, ra=ra, dec=dec, dir_cos=dir_cos,$
;         x_ra_rad=x_ra_rad ,x_dec_rad=x_dec_rad ,bad_data=bad_data
;
; CALLED BY:
;
;
; CALLS:
;	none
;
; INPUTS:
;
;	ut - time in seconds from 79/1/1 corresponding to vectors x_ra, x_dec
;	x_ra - GRO X-axis pointing direction right ascension in .0001 RAD
;	x_dec -GRO X-axis pointing direction declination in .0001 RAD
;	z_ra - GRO Z-axis pointing direction right ascension in .0001 RAD
;	z_dec -GRO Z-axis pointing direction declination in .0001 RAD
;	pnt_ndx - pointing index to start of new pointing direction, ref. to ut
;
; OPTIONAL INPUTS:
;	none
;
; OUTPUTS:
;       none explicit, only through commons;
;
; OPTIONAL OUTPUTS:
;	cos_det - direction cosines for each of BATSE LADs for each time sample
;	   cos_det is dimensioned fltarr(8,n_elements(ut))
;	cos_spec - direction cosines for each of BATSE SPECs for each time sample
;	   cos_spec is dimensioned fltarr(8,n_elements(ut))
;       bad_data - set to 1 if data quality is too poor
; COMMON BLOCKS:
;	none
;
; SIDE EFFECTS:
;	none
;
; RESTRICTIONS:
;	none
;
; PROCEDURE:
;	none
;
; MODIFICATION HISTORY:
; 	written, ras, sometime 1994
;	mod ras, 22-jun-1995, crossp=>f_crossp
;	mod ras, 8-sep-1995, added cos_spec
;
;-
;******************************
pro det_point, ut, x_ra, x_dec, pnt_ndx, cos_det,z_ra,z_dec,$
	celestial_fov=batsefov, cos_spec=cos_spec, $
	spacecraftfov=batse,x2fov=x2fov,ra=ra,dec=dec,dir_cos=dir_cos,$
	x_ra_rad=x_ra_rad,x_dec_rad=x_dec_rad,bad_data=bad_data

bad_data=0 ;it is assumed that at least one valid pointing direction is passed
	   ; bad_data=1 means that only bad data has been passed
;Convert X_RA and X_DEC from .0001 rad stored in integer*2 format
;to floating point in radians
;X_RA goes from 0 to 10000*2*!pi, unsigned
;X_DEC goes from -10000*!pi/2 to 10000*!pi/2 signed
	nx=n_elements(x_ra)
	;Interpret x_ra as the low 2 bytes of a longword
	x_ra_rad= (x_ra and 65535l)/1e4 ;ra in radians
	x_dec_rad=x_dec/1e4 ;declination in radians
	xc=transpose(sphcart(x_ra_rad*!radeg,x_dec_rad*!radeg)) ;3 by nx
	z_ra_rad= (z_ra and 65535l)/1e4 ;ra in radians
	z_dec_rad=z_dec/1e4 ;declination in radians
	zc=transpose(sphcart(z_ra_rad*!radeg,z_dec_rad*!radeg)) ;3 by nx
	yc=f_crossp(zc,xc)	
;***************************
;GET THE POSITION OF THE SUN IN CELESTIAL CARTESIAN COORDINATES
;Get Ra and Dec of the Sun at center of orbit
	sradec=solephut(avg(ut(where(ut gt utime('90/01/01')))))
;Vector to Sun in cartesian coordinates
	sunxyz=sphcart(sradec(0),sradec(1)) 
	;sunxyz is (1,3)
;******************************
;LOOK FOR BAD SC RA AND DEC - RA=0 AND DEC=0
;			 - RA > 2*!PI, ABS(DEC) > !PI/2
;			 - SUN > 42 DEGREES FROM THE XZ PLANE 
;			 - X AND Z NOT PERPENDICULAR
	wbad=where( ( (x_ra eq 0) and (x_dec eq 0)) or $
		    ( (z_ra eq 0) and (z_dec eq 0)) or $
		    ( x_ra_rad gt 1.01*2*!pi) or $
		    ( z_ra_rad gt 1.01*2*!pi) or $
		    ( abs(x_dec_rad) gt 1.01*!pi/2 ) or $
		    ( abs(z_dec_rad) gt 1.01*!pi/2 ) or $
$;		     X AND Z NOT PERPENDICULAR
		    ( abs( [1.,1.,1.]#(zc*xc) ) gt .01 ) or $ ;x perp. z
$; 		     SUN > 42 DEGREES FROM THE XZ PLANE 
		    ( abs(sunxyz#yc) gt cos((90.-42.-2.)/!radeg)  ) )
if !err eq nx then begin
	bad_data=1
	return
endif
if !err ge 1 then begin
;FILL THE BAD INDICES WITH THE PRECEDING GOOD INDEX
	;FIND THE GROUPS OF BAD INDICES
	group_intervals,indices=wbad,groups=groups ;groups is 2 by ng 
	; ng is the number of separate groups of consecutive samples in wbad
		ng=n_elements(groups)/2
		for i=0,ng-1 do begin 
	;fill bad pointing with previous good pointing
		  good=groups(0,i)-1
			if good lt 0 then good=groups(1,i)+1 ;take it after
		  x_ra_rad(groups(0,i))= $
		   x_ra_rad(good)+fltarr(groups(1,i)-groups(0,i)+1)
		  x_dec_rad(groups(0,i))= $
		   x_dec_rad(good)+fltarr(groups(1,i)-groups(0,i)+1)
		endfor
endif ;end fill for bad spacecraft pointing
	pnt_ndx=[0]
	istart=0 ;calculate pointing at start of time sample
	new_pointing: ;LOOK FOR CHANGES IN POINTING 
		any_change,istart,x_ra_rad,x_dec_rad,ichange
		 if ichange gt istart then begin
			 pnt_ndx=[pnt_ndx,ichange]
			 istart=ichange
			 if istart lt n_elements(x_ra)-2 then goto,new_pointing
		endif 
		nchange=n_elements(pnt_ndx)
;;;;;;;;;;;;;;;;;;;;;;; Detectors' Aspects in Spacecraft Coordinates
;given the spacecraft z-axis as straight up along COMPTEL and the x-axis
;is pointed in the OSSE scan plane, find the direction cosines of the normals
;to the 8 BATSE LAD's
	batse=fltarr(3,8)
	index=indgen(1,8)
	for j=0,2 do batse(j,0)=sqrt(1/3.)*(-1)^(index/2^(2-j))
;;;;;;;;;;;;;;;;;;;;;;;;;
;ADD SPECTROSCOPY DETECTORS POINTING
	zcos=batse(2,*)
	zang = acos(zcos)*!radeg
	zang2=zang+(-1)^index*19
	zsin= sin(zang/!radeg)
	zcos2 = cos(zang2/!radeg)
	zsin2 = sin(zang2/!radeg)
	lad2spec = [ zsin2/zsin, zsin2/zsin, zcos2/zcos]
	spec = batse*lad2spec
;;;;;;;;;;;;;;;;;;;;;;;;;
;CREATE TRANSFORMATIONS FROM DETECTOR TO CELESTIAL COORDINATES
	batsefov=fltarr(3,8,nchange) ;1 matrix for every pointing change
	specfov =fltarr(3,8,nchange) ;1 matrix for every pointing change
	for j=0,nchange-1 do begin
	  pj=pnt_ndx(j) 
	  x2fov=(reform([xc(*,pj),yc(*,pj),zc(*,pj)],3,3))
	  fov2x=transpose(x2fov)
	  ;detector directions in celestial
	  batsefov(0,0,j)=x2fov#batse
	  specfov(0,0,j)=x2fov#spec
	endfor
; (vector in celestial coord.)=(vec. in detector coor)# x2fov
;****************
;COMPUTE THE DIRECTION COSINES TO THE SUN FOR EACH LAD FOR EACH SAMPLE 
	nut=n_elements(ut)
	npt=[pnt_ndx,nut]
	cos_det=fltarr(8,nut)
	cos_spec=fltarr(8,nut)
	for j=0,nchange-1 do begin
	  dir_cos=reform(sunxyz#reform(batsefov(*,*,j))) ;tran(1x3 # 3x8)=8
	  cos_det(0,npt(j))=dir_cos#replicate(1,1,npt(j+1)-npt(j))
	  dir_cos2=reform(sunxyz#reform(specfov(*,*,j))) ;tran(1x3 # 3x8)=8
	  cos_spec(0,npt(j))=dir_cos2#replicate(1,1,npt(j+1)-npt(j))
	endfor
return &end