pro WPHPLANE, TFtype,basis,pkt,titlestr=titlestr,nTFR=nTFR,qmf=qmf
;+
; NAME:
;	WPHPLANE
;
; PURPOSE:
;	This procedure partitions phase space by rectangular blocks. 
;
; CATEGORY:
;	Wavelets.
;
; CALLING SEQUENCE:
;	WPHPLANE, TFtype,btree,pkt,titlestr=titlestr,nTFR=nTFR,qmf=qmf 
;
; INPUTS:
;    TFtype =    string: type of TF-packets used ('WP','CP') 
;    btree  =    basis tree 
;    pkt    =    wavelet or cosine packet table 
;    titlestr  = signal name 
;    nTFR   =    number of x-points in phaseplane image (optional)
;    qmf    =    qmf for calculating WP phase plane location (optional)
;
; OUTPUTS:
;    an image plot with colored rectangles based on  
;    recursive dyadic partitioning of y axis according 
;    to splits in basis tree 
;
; EXAMPLE:
;  To create a ``Heisenberg image" plot of data analyzed with WPANALYS,
;  make the following call:
;
;	IDL>WPHPLANE,`WP',btree,wp,titlestr = st, nTFR = 64
;
;  where `WP' is wavelet packet, btree and wp are outputs from WPANALYS,
;  st is the string identifier for the data (`seismic', `sunspots', etc.),
;  nTFR is the number of boxes (dyads!) on a side. Keep nTFR on the low side, 
;  If nTFR = 256, or higher you'll be waiting 5 minutes or more for output.
;
; NOTE: 
;	This routine also performs a Cosine-Packet Analysis ('CP'), however,
;	cosine packet analysis isn't currently a function fully implemented
;	in Wavelet Workbench. So consider this 'CP' option a 'hook' for now.
;
;
; MODIFICATION HISTORY:
; 	Written by:	Amara Graps		November, 1995
;Translated from MatLab Wavelab routine: imagephaseplane.m
;-

;Color of text strings
IF !p.background EQ 0 THEN cl = 255 ELSE cl = 0

;Set parameters if keywords are not set

IF (KEYWORD_SET(titlestr) EQ 0) THEN titlestr = ' '
IF (KEYWORD_SET(nTFR) EQ 0) THEN nTFR = 128
IF (KEYWORD_SET(qmf) EQ 0) THEN Corrected = 0 ELSE Corrected = 1
;Note Corrected = 0 = False and Corrected = 1 = True

t = SIZE(pkt) 
L = t(1)		;number of columns
n = t(2)		;number of rows
 
; initialize tree traversal stack  
dstack = FLTARR(2^L) 
bstack = FLTARR(2^L) 
k = 1
dstack(k-1) = 0 
bstack(k-1) = 0 

ss = WNNORM(pkt(0, *)) 
 
TFPLane = FLTARR(nTFR,nTFR) 

WHILE (k GT 0) DO BEGIN 
	d = dstack(k-1) 
	b = bstack(k-1) 
 	k=k-1 
 	
	IF ((basis(WNODE(d,b)-1) EQ 0) ) THEN BEGIN  
		; terwmminal node 
		ylo = b/2^d 
 		yhi = (b+1)/2^d 
		nylo = 1 + FLOOR(nTFR * ylo) 
		nyhi = 1 + FLOOR(nTFR * yhi) 
		
		;use the IDL way to create an index array
		aa = nylo 	;val of 1st element
		bb = 1	;increment
		cc = nyhi	;val of last element
		indx_array = WMKIARRY( aa, bb, cc)+1		

		yind = WPAROUND(indx_array,nTFR);
		coeffs = n * (pkt(d,WPACKET(d,b,n)) /ss) ^2 
		nbox = n / 2^d   ; = length(coeffs) 

		Checktype =  STRCOMPRESS(TFtype,/remove_all)	;remove white space
		Checktype = STRUPCASE(Checktype)	;make upper case

		CASE 1 OF		;1 = True 
			((Checktype EQ 'WP') AND (Corrected)): BEGIN
				pos = WWPLOC(d,b,0,qmf,n)
				END
			ELSE: pos = 0
		ENDCASE

				
		FOR  j=0, (nbox-1) DO BEGIN 
			jnew = WREM((pos/n)*nbox + j, nbox) 
			
			IF ( coeffs(j) GT 0) THEN BEGIN 
				nxlo = 1 + FLOOR((nTFR *jnew)/nbox) 
				nxhi = 1 + FLOOR((nTFR *(jnew+1))/nbox) 
				
				;use the IDL way to create an index array
				aa = nxlo 	;val of 1st element
				bb = 1	;increment
				cc = nxhi	;val of last element
				indx_array = WMKIARRY( aa, bb, cc)+1	

				xind = WPAROUND(indx_array,nTFR);
				;Note, the "-1" is handled in WMKIARRY, so it's not necessary in TFPlane

				CASE 1 OF		;1 = True 
					(Checktype EQ 'CP'): BEGIN
					 	;Note: Matlab does the following line (IDL can't do this)
						;TFPlane(yind-1, xind-1) = TFPlane(yind-1, xind-1) + coeffs(j) 
						;Assign xind and yind to TFPlane
						FOR i = 0, N_ELEMENTS(xind)-1 DO BEGIN
							TFPlane(yind-1,xind(i)-1) = TFPlane(yind-1,xind(i)-1) + coeffs(j) 
							END

						END
					ELSE: BEGIN
						;Note: Matlab does the following line (IDL can't do this)
						;TFPlane(xind-1, yind-1) = TFPlane(xind-1, yind-1) + coeffs(j) 
						;Assign xind and yind to TFPlane
						FOR i = 0, N_ELEMENTS(yind)-1 DO BEGIN
							TFPlane(xind-1,yind(i)-1) = TFPlane(xind-1,yind(i)-1) + coeffs(j) 
							END


						END
				ENDCASE

			END ;IF
			
		END ;FOR  j

		
		ENDIF ELSE BEGIN 		; endif ((basis(WNODE(d,b)-1) eq 0) ) 
			k = k+1 
			dstack(k-1) = d+1 
			bstack(k-1) = 2*b 
			k = k+1 
			dstack(k-1) = d+1 
 			bstack(k-1) = 2*b+1 
		END 	;ENDELSE
		
	END ;WHILE (k GT 0)
;  

LOADCT, 5	
ptitle = 'Phase plane: ' + titlestr

TFMax = MAX(MAX(TFPlane)) 
TFPlane = TFPlane * ( 256/TFMax) 

;colormap(1-gray(256)) ;Matlab WaveLab does this, it is necessary for IDL?

WIMAGE,TFPLane, xrange= WLNSPACE(0,1,nTFR), yrange=WLNSPACE(0,1,nTFR), $
	title=ptitle, xtitle='Time', ytitle='Frequency' 

;Note I'm not using the following Matlab commands (don't think that they
;are necessary, but not 100% sure. --AG).
;axis('xy') ;causes the graph to go to cartesian coord form 
;bottom-to-top y numbering, left-to-right x numbering
;axis ([0 1 0 1 ]) ;sets the scaling from 0 to 1 (x and y) on 
;current plot
 
; 
; Algorithm Source:  WaveLab Version 0.600
; WaveLab WWW site: http://playfair.stanford.edu/ 
; WaveLab Questions? e-mail wavelab@playfair.stanford.edu
;    
     
 
END  
;***************************************************