;+
; NAME:
;	FRACTAL_CLOUD
; PURPOSE:
;
;	Generate a fractal cloud as a density distribution map,
;	returning a 3D array (1D or 2D if requested). The density map
;	is generated by binning points generated by a randomn recursive
;	fractal algorithm, as described by Elmegreen, 1997, ApJ 477: 196-203.
;
; CALLING:
;
;	dens_map = Fractal_Cloud( NDIM=, NLEVELS=, NP_FACTOR=, LSCALE=, $
;					SIZE_MAP=, /VERBOSE, /DISPLAY )
;
; KEYWORD INPUTS:
;
;	NDIM = dimension of space in which to embed fractal, default is 3.
;
;	NLEVELS = number of LEVELS in fractal heiarchy, no default, must give.
;
;	NP_FACTOR = integer factor (>1) setting the number of random points
;		generated at each new level around and replacing each single
;		point of previous level, default is 12.
;		Total number of random points generated = NP_FACTOR^NLEVELS.
;
;	LSCALE = the spatial scale divisor from one level to next.
;		The spatial extent of new NP_FACTOR random points at next level
;		is decreased by divisor LSCALE from previous level,
;		thus spatial scale at level N = LSCALE^(-N).
;		Can be scalar or vector with NDIM elements, except
;		in case of spherical symmetry only scalar is used. Default=3.
;
;	SIZE_MAP = dimensions (resolution in pixels) of output density map,
;		default is [128,128,128].
;
;	/SPHERICAL : generate points with spherical symmetry instead of cubic.
;
;	RADIUS_EXP = exponent applied to random radii, for /SPHERICAL case only,
;		default = 1/3 (for uniform in 3D sphere).
;
;	/SAVE_RS : save the currently generated random number sequences into
;		IDL/XDR save files (in current directory) which can then be
;		later restored and used to obtain the same fractal realization
;		with different dimension (by changing LSCALE) or resolution.
;
;	DIR_SAVED_RS = optional string specifying directory containing IDL/XDR 
;	    save files in which previously generated random numbers are stored.
;		For UNIX, must use a "/" character at then end of dir string.
;
;	RANGE_XYZ = float array of dimensions (2,NDIM) specifying min-max
;		range of each dimension for binning the points. Default range
;		is a sphere which contains all points.
;		
;	/VERBOSE
;	/DISPLAY
;
; KEYWORD OUTPUTS:
;
;	RANGE_XYZ = min-max range used for each dimension in binning the points.
;
;	DHISTOGRAM = histogram of the density map.
;	AVG_DENS = average density of the cloud pixels/voxels, excluding ICM.
;
; OUTPUT:
;	Function returns density distribution map,
;	an integer array of the requested dimension and size.
;
; EXTERNAL CALLS:
;	function Bin3D
;	function Bin2D
;	function Histo
;	function VarType
;	function trapez
;	pro get_window
;	pro tvs
; PROCEDURE:
;	Use the uniform random number generator to make points at each Level,
;	which are rescaled in a descending heirarchy, a modified version of
;	the algorithm described by Elmegreen, 1997, ApJ 477: 196-203.
;	The final NP_FACTOR^NLEVEL points are binned into a density map array,
;	but not all at once, just one group at a time, looping over NP_FACTOR.
; HISTORY:
;	Written: Frank Varosi, HSTX @ NASA/GSFC, 1997.
;-

function Fractal_Cloud, NDIM=Ndim, NLEVELS=NLev, NP_FACTOR=Npf, LSCALE=Lscale, $
			DIR_SAVED_RS=dirseq, SPHERICAL=spheric, SIZE_MAP=nvox, $
			VERBOSE=verbose, DISPLAY=display, LOG=minLog, $
			DHISTOGRAM=dmh, AVG_DENS=avg_dens, RANGE_XYZ=rxyz, $
			RADIUS_EXP=rexp, SAVE_RS=save_rs

if N_elements( NLev ) ne 1 then begin
	print,"syntax:	dens_map = Fractal_Cloud( NLEVEL=, NP_FACTOR=, LSCAL=,$"
	print,"				NDIM=, SIZE_MAP=, /VERBOSE, /DISPLAY )"
	return,0
   endif

if N_elements( Npf ) ne 1 then Npf = 12
if N_elements( Ndim ) ne 1 then Ndim = 3

if keyword_set( spheric ) then begin
	if N_elements( Lscale ) gt 1 then Lscale = Lscale(0)
	if N_elements( Lscale ) ne 1 then Lscale = 3
	if N_elements( rexp ) ne 1 then rexp = 1./3.
	Ndim = 3
 endif else begin
	if N_elements( Lscale ) eq 1 then Lscale = replicate( Lscale, 3 )
	if N_elements( Lscale ) LE 0 then Lscale = replicate( 3.0, 3 )
  endelse

Lscale = float( Lscale )
slice = 0.01

if VarType( dirseq ) eq "STRING" then begin

	restore, VERB=keyword_set( verbose ), dirseq + "RanSeq_L=1_k=00.idl"
	szr = size( ranseq )
	if szr(1) LT szr(2) then ranseq = transpose( ranseq )
	szr = size( ranseq )
	Npf = szr(1)
	Ndim = szr(2)
	if keyword_set( verbose ) then help, ranseq, Level
	restrseq = 1

	if keyword_set( spheric ) and (szr(2) ne 3) then begin
		message,"need 3D ran.seq. for spherical option",/INFO
		help, ranseq
		return,0
	   endif

 endif else begin

	ranseq = randomu( seed, Npf, Ndim )

	if keyword_set( save_rs ) then begin
		Level = 1
		save, ranseq, Level, FILE="RanSeq_L=1_k=00.idl",/VERB
	   endif
  endelse

if keyword_set( spheric ) then begin

	w = 2*ranseq(*,0)-1
	c = sqrt( 1 - w*w )
	a = !PI * (2*ranseq(*,1)-1)
	xo = [ [ c * cos( a ) ],[ c * sin( a ) ],[ w ] ]
	w=0 & c=0 & a=0
	radii = ( ranseq(*,2)^rexp )/Lscale
	for i=0,Ndim-1 do xo(*,i) = xo(*,i) * radii

 endif else begin

	xo = 2*( ranseq - 0.5 )
	Lscale = Lscale(0:Ndim-1)
	for i=0,Ndim-1 do xo(*,i) = xo(*,i)/Lscale(i)
  endelse

if keyword_set( display ) then get_window,0,/show
D_cloud = aLog( Npf )/( total( aLog( Lscale ) )/N_elements( Lscale ) )

if keyword_set( verbose ) then begin
	help, Npf, NLev, Ndim, rexp, D_cloud
	print," Lscale:",Lscale
   endif

for Lev = 2, NLev-1 do begin		;Loop over heirarchy of Levels:

	szx = size( xo )
	Nx = szx(1)
	xn = fltarr( Nx, Npf, Ndim )

	for k=0,Npf-1 do begin

		if keyword_set( restrseq ) then begin

			restore, dirseq + "RanSeq_L=" + strtrim( Lev, 2 ) $
			    + "_k=" + strmid( strtrim( k+100, 2 ),1,2 ) + ".idl"
			szr = size( ranseq )
			if szr(1) LT szr(2) then ranseq = transpose( ranseq )

		 endif else begin

			ranseq = randomu( seed, Nx, Ndim )

			if keyword_set( save_rs ) then begin
				Level = Lev
				file = "RanSeq_L=" + strtrim( Level,2 )+"_k="+ $
					strmid( strtrim( k+100, 2 ),1,2 )+".idl"
				save, ranseq, Level, FILE=file
			   endif
		  endelse

		if keyword_set( spheric ) then begin

			w = 2*ranseq(*,0)-1
			c = sqrt( 1 - w*w )
			a = !PI * (2*ranseq(*,1)-1)
			x = [ [ c * cos( a ) ],[ c * sin( a ) ],[ w ] ]
			w=0 & c=0 & a=0
			radii = ( ranseq(*,2)^rexp ) * Lscale^(-Lev)
			for i=0,2 do x(*,i) = x(*,i) * radii
			radii = 0
			xn(*,k,*) = xo + x
			x = 0

		 endif else for i=0,Ndim-1 do xn(*,k,i) = xo(*,i) + $
					2*(ranseq(*,i) - 0.5) * Lscale(i)^(-Lev)
		ranseq=0
	  endfor

	xo = 0
	xo = reform( xn, Nx*Npf, Ndim )
	xn = 0
	if keyword_set( verbose ) then help,Lev,xo

	if keyword_set( display ) then begin
	   CASE Ndim OF
		1:	plot, hv, histo( xo, hv ), PS=10
		2:	plot, xo(*,0), xo(*,1), ps=3, XR=[-1,1], YR=[-1,1]
		3: BEGIN
			w = where( abs( xo(*,2) ) LE slice, nw )
			if (nw GT 1) then plot, xo(w,0), xo(w,1), ps=3, $
						XR=[-1,1], YR=[-1,1]
		     END
	     ENDCASE
		empty
	   endif
  endfor				;end Loop over Levels (except Last one).

radxo = max( vec_norm( reform( xo ) ) ) + min(Lscale)^(-Lev)
sz = size( rxyz )

if sz(0) ne 2 then begin

	rxyz = fltarr( 2, Ndim )

	for i=0,Ndim-1 do begin
		rxyz(0,i) = min( xo(*,i), MAX=m )
		rxyz(1,i) = m
	  endfor

	if keyword_set( verbose ) then begin
		print," range_xyz: "
		print, rxyz
	   endif

	for i=0,Ndim-1 do rxyz(*,i) = [-radxo,radxo]
   endif

if keyword_set( verbose ) then begin
	print," using range_xyz: "
	print, rxyz
	if sz(0) eq 2 then help, radxo
   endif

if N_elements( nvox ) LE 0 then nvox = 128
if keyword_set( verbose ) then print," Nvox: ",nvox
xn = xo

for k=0,Npf-1 do begin		;Bin the final Level as points are generated:

	if keyword_set( restrseq ) then begin

		restore, dirseq + "RanSeq_L=" + strtrim( NLev, 2 ) $
			+ "_k=" + strmid( strtrim( k+100, 2 ),1,2 ) + ".idl"
		szr = size( ranseq )
		if szr(1) LT szr(2) then ranseq = transpose( ranseq )

	 endif else begin 

		szx = size( xo )
		ranseq = randomu( seed, szx(1), szx(2) )

		if keyword_set( save_rs ) then begin
			Level = NLev
			file = "RanSeq_L=" + strtrim( Level,2 ) + "_k=" + $
				strmid( strtrim( k+100, 2 ),1,2 ) + ".idl"
			save, ranseq, Level, FILE=file
		   endif
	  endelse

	if keyword_set( spheric ) then begin

		w = 2*ranseq(*,0)-1
		c = sqrt( 1 - w*w )
		xn(*,2) = w
		w = 0
		a = !PI * (2*ranseq(*,1)-1)
		xn(*,0) = c * cos( a )
		xn(*,1) = c * sin( a )
		c=0 & a=0
		radii = ( ranseq(*,2)^rexp ) * Lscale^(-NLev)
		for i=0,2 do xn(*,i) = xn(*,i) * radii
		radii = 0
		xn = xo + temporary( xn )

	 endif else for i=0,Ndim-1 do xn(*,i) = xo(*,i) + $
					2*( ranseq(*,i)-0.5 )* Lscale(i)^(-NLev)
	ranseq=0

	CASE Ndim OF

	 3:	dens_map = Bin3D( xn, NVOX=nvox, TYPE=2, VOXEL_DENS=dens_map, $
				XRAN=rxyz(*,0), YRAN=rxyz(*,1), ZRAN=rxyz(*,2) )

	 2:	dens_map = imagxy( xn, NPIX=nvox, TYPE=2, IMAGE_DENS=dens_map, $
						XRAN=rxyz(*,0), YRAN=rxyz(*,1) )

	 1: BEGIN
		if( k eq 0 ) then dens_map = Histo( xn, NBIN=nvox, $
						MIN=rxyz(0), MAX=rxyz(1) ) $
			else dens_map = dens_map + Histo( xn, NBIN=nvox, $
						MIN=rxyz(0), MAX=rxyz(1) )
	      END

	 else:	message,"NDIM = " + strtrim( Ndim, 2 ) + " not supported",/INFO
	 ENDCASE

	if (k eq 0) and keyword_set( verbose ) then help,dens_map

	if keyword_set( display ) or $
	   keyword_set( verbose ) or (k eq (Npf-1)) then begin
		dmh = histogram( dens_map )
		ff_icm = dmh(0)/total( dmh )
		max_dens = N_elements( dmh )-1
		hv = 1 + findgen( max_dens )
		hd = dmh(1:*)
		avg_dens = trapez( hd * hv, hv )/trapez( hd, hv )
		if keyword_set( verbose ) then $
		 print, k, max_dens,"=max", avg_dens,"=avg", ff_icm,"=ff_icm", $
			radxo - max( vec_norm( xn ) )
	   endif

	if keyword_set( display ) then begin

	    CASE Ndim OF

	     3: BEGIN

		for i=1,Ndim do begin
			CASE i OF
			  1:	image = reform( dens_map(nvox(0)/2,*,*) )
			  2:	image = reform( dens_map(*,nvox(1)/2,*) )
			  3:	image = dens_map(*,*,nvox(2)/2)
			 ENDCASE
			get_window,i,/SHOW, XS=340, YS=256, $
					XP=!DEVX-700, YP=!DEVY-290*i+55
			tvs, image, LOG=minLog, COL=2,/ERAS
		  endfor

		for i=1,Ndim do begin
			get_window,10+i,/SHOW, XS=350, YS=256, $
					XP=!DEVX-350, YP=!DEVY-290*i+55
			tvs, total( dens_map, i ), COL=2, /ERAS, LOG=minLog
		  endfor
	      END

	     2: BEGIN
			get_window,1,/SHOW, XS=640, YS=512
			tvs, dens_map, COL=2, /ERAS, LOG=minLog
	          END

	     1: BEGIN
			get_window,1,/SHOW
			plot, dens_map, PS=10
	          END

	     else:
	    ENDCASE

		get_window,0,/SHOW, XS=!DEVX-700, YS=!DEVX-700, $
				XP=0, YP=!DEVY-!DEVX+700
		plot, hv, hd > 0.1, /YTYPE, PS=10
		oplot, hv, exp( -hv/avg_dens ) * hd(0) * exp(1/avg_dens),LINE=2
		empty
	   endif
  endfor

if keyword_set( verbose ) then begin
	print," # points =", total( dens_map ), Long( Npf )^NLev
	help, Npf, NLev, rexp
	print," Lscale =",Lscale
	cdmax = float( nvox(0) )^3
	help, ff_icm, 1 - cdmax^((D_cloud/3)-1), D_cloud
   endif

return, dens_map
end