;+
; NAME:
;	catimg
; PURPOSE:
;	To render 3D images of CATS density matrices.
; CALLING SEQUENCE:
;	catimg,AX,AZ,[RES=,THRESH=,SCALE=,CUT_THRESH=,SLCTYP=,PER=,TYPE=,/VIDEO]
; INPUTS:
;	AX	Rotate image about X axis in degrees.
;	AZ	Rotate image about Z axis in degrees.
; KEYWORDS:
;
;	RES	A scalar specifying the resolution of image in square pixels.
;	THRESH	Contouring threshold.
;	SCALE	Scaling factor for size of volume relative to coordinates
;	CUT_TRESH	> 0 then object is cut in half and then rendered.
;	SLCTYP	If CUT_THRESH > 0. then this determines what kind of slice
;		to make. 
;		+ or - 1 = Positive or Negative X-Axis slice.
;		+ or - 2 = Positive or Negative Y-Axis slice.
;		+ or - 3 = Positive or Negative Z-Axis slice.
;	PER	Scalar which determines perspective distance.
;	VIDEO	if set then object is rendered with a black background and
;		white double-thick axes and labelling.
;	TYPE	= 0 40x10x10 matrix where the unit cell is (.03 AU/cell)^3
;		= 1 40x20x10 matrix where the unit cell is (.03 AU/cell)^3
;		= 2 40x20x10 matrix where the unit cell is (.05 AU/cell)^3
;		    (For May 24th event only)
; OUTPUTS:
; COMMON BLOCKS:
;	CATDATA
; RESTRICTIONS:
; PROCEDURE:
;	Use CATIMG to produce TIFF images.
; EXAMPLES:
;	CATIMG,30,160,thresh=60.,res=300,/VIDEO,file='may7_50.tif',type=0
;
; This sequence reads the may7_50.dat CATS density file and produces a 3D
; contour image with a threshold value 60, 30 deg. above X axis,
; 160 degrees around Z axis in a 300x300 pixel image in the video format in a 
; TIFF file called 'may7_50.tif'.
; 
; MODIFICATION HISTORY:
; - Tom Davidson May 1993 (UCSD)
;-
;scale = [7.0177130E-02,0.1037752,0.1465857,0.1938648,0.2411200,0.2883753, $
;        0.3356305,0.3828857,0.4301649,0.4776014,0.5250890,0.5725765,0.6200640,$
;        0.6675754,0.7150630,0.7625505,0.8100380,0.8575255,0.9050370,0.9525245,$
;        1.000012,1.047500,1.094987,1.142497,1.189983,1.237470,1.284956,$
;	 1.332442,1.379952,1.427438,1.474925,1.522411,1.569898,1.617408,$
;	 1.664894,1.712380,1.759867,1.807353,1.854863,1.902349,1.949836]

;scale = [0.169,0.221,0.292,0.381,0.470,0.558,0.647,0.735,0.824,0.912, $
;         1.00,1.088,1.177,1.265,1.354,1.443,1.533,1.622,1.712,1.802,1.891]

pro CATIMG,INFILE,ax=AX,az=AZ,	$
	res=RES,thresh=THRESH,scale=SCALE,cut_thresh=CUT_THRESH,	$
	slctyp=SLCTYP,per=PER,video=VIDEO,outfile=OUTFILE,	$
	type=TYPE,dev=DEV,anim=ANIM, prompt=PROMPT

ANIM   = keyword_set(ANIM)
PROMPT = keyword_set(PROMPT)
PROMPT = PROMPT and not ANIM

NX = n_elements(AX)
NZ = n_elements(AZ)
if NX gt 1 and NZ gt 1 then begin
    if NX ne NZ then message, 'AX and AZ must have same size'
endif else if NX gt 1 then begin
    if NZ eq 0 then AZ = 0
    AZ = replicate(AZ,NX)
endif else if NZ gt 1 then begin
    if NZ eq 0 then AX = 0
    AX = replicate(AZ,NZ)
endif else if NX+NZ eq 1 then begin    
    if n_elements(AX) eq 0 then AX = 30
    if n_elements(AZ) eq 0 then AZ = 30
endif else if NX+NZ eq 0 then begin
    if ANIM then begin
	AX = replicate(10,36)
	AZ = -135+10*indgen(36)
    endif else begin
	AX = 30
	AZ = 30
    endelse
endif
NX = n_elements(AX)

if n_elements(SCALE) eq 0 then SCALE = 1.
NS = n_elements(SCALE)

VIDEO = keyword_set(VIDEO)
if n_elements(INFILE) eq 0 then INFILE = filepath(root=getenv('DAT'),'may24i.dat')
if n_elements(THRESH) eq 0	then THRESH = 50.
if n_elements(CUT_THRESH) eq 0	then CUT_THRESH = 0.
if n_elements(SLCTYP) eq 0	then SLCTYP = 1
if n_elements(RES) eq 0		then RES = 400
if n_elements(OUTFILE) eq 0	then OUTFILE = 'test.tiff'
if n_elements(TYPE) eq 0	then TYPE = 2
if n_elements(DEV) eq 0		then DEV = 'X'

if PROMPT then begin
    echo,'Input file name',INFILE
    echo,'Angle of rotation about X axis:',AX
    echo,'Angle of rotation about Z axis:',AZ
    echo,'Resolution of image in pixels:',RES
    echo,'Isosurface threshold value:',THRESH
    echo,'Axis type:
    echo,'   0=40x10x10 (May 7),'
    echo,'   1=40x20x10 (May 7),'
    echo,'   2=40x20x10 (May 24) ):',TYPE

    echo,'Image type:
    echo,'   0=Black on White,'
    echo,'   1=White on Black:',VIDEO
    echo,'Plot Type (X,PCL,TIFF)',DEV
    DEV = strupcase(DEV)
endif

; Read file into dummy array.
; The CAT program stores the array in a rather screwy way.
; If the density array has dimension SzX,SzY,SzZ then bin_read
; returns a 2D array DUM(SzY*SzX,SzX).
; SzZ is specified to reform the array to its proper size.

SzZ = 10

if not bin_read(INFILE,DUM,/bor,/afltarr,error=S) then message, S

S = size(DUM)  &  SzY = S(1)/SzZ  &  SzX = S(2)
DUM = reform(DUM,SzZ,SzY,SzX)
SDEN = fltarr(SzX,SzY,SzZ,/nozero)
for i=0,SzX-1 do for j=0,SzY-1 do SDEN(i,j,*) = DUM(*,j,i)

DUM = 0.
SDEN( where(SDEN lt 1.0) ) = .0

whatis, SDEN
message, 'Total in SDen ='+string(total(SDen)), /info

XOff = 10
YOff =  0
DIMX = 20
DIMZ = 10
case TYPE of
0: DIMY = 10			; 40x10x10
1: DIMY = 20			; 40x20x10
2: DIMY = 20			; 40x20x10 (May 24 CME)
else: message, 'Type is invalid'
endcase

VOL = SDEN(XOff:XOff+DIMX-1,YOff:YOff+DIMY-1,0:DIMZ-1)

;  For best results you may have to customize these in order to get
;  the best possible half-shells.

if CUT_THRESH gt 0. then begin
   case SLCTYP of
    1:	vol=vol(12:19, 0:9, 0:9 )
   -1:	vol=vol( 0: 9, 0:9, 0:9 )
    2:	vol=vol( 0:19, 5:9, 0:9 )
   -2:	vol=vol( 0:19, 0:4, 0:9 )
    3:	vol=vol( 0:19, 0:9, 5:9 )
   -3:	vol=vol( 0:19, 0:9, 0:4 )
   else:message, 'SLCTYP is invalid.'
   endcase
endif

for I=0,NX-1 do for J=0,NS-1 do		$
	catrend, VOL,THRESH,SCALE(J),CUT_THRESH,SLCTYP,AX(I),AZ(I),1,	$
		RES,PER,VIDEO,TYPE,OUTFILE,DEV

return  &  end