;+
; NAME:
;	gridgen
; PURPOSE:
;	Generate array of x,y,z coordinates for a multi-dim area
; CALLING SEQUENCE:
	FUNCTION gridgen, nn, yy, zz, $
		unit	= unit		, $
		merge	= merge		, $
		origin	= origin	, $
		edge	= edge		, $
		open	= open		, $
		one 	= one		, $
		range	= range 	, $
		multid	= multid	, $
		double	= double
; INPUTS:
;	n			array[ndim]; type: integer
;				(ndim = 1,2,3)
;				dimensions of area in x, x-y, or x-y-z directions
; OPTIONAL INPUT PARAMETERS:
;	/one		if set, scales output array to range [0,1] in all dimensions
;	one=one		array[ndim]; type: integer with value 0 or 1
;				set scaling for each of the dimension in 'n' separately
;
;	/open		if set, then the pixel centers for the 'pixel edges' [0,n]
;				is returned
;	open=open	array[ndim]; type: integer with value 0 or 1
;				set 'open' status for each of the dimensions separately
;				This is the same as gridgen1d(n, origin=-0.5)
;				The origin=0.5 is added to the origin specified in keyword
;				origin.
;
;	/edge		if set, then coordinates for the 'pixel edges' of pixels
;				[0,..n-1] is returned
;	edge=edge	array[ndim]; type: integer with value 0 or 1
;				set 'edge' status for each of the dimensions separately
;				This is the same as gridgen(n+1, origin=0.5)
;				The origin=0.5 is added to the origin specified in keyword
;				origin.
;
;				Example: gridgen(5, /edge) = [-0.5,0.5,1.5,2.5,3.5, 4.5]
;	origin=origin
;				array[ndim]; type: float
;				defines the origin in units of output grid
;				(i.e. usually in array indices; but if /one or range are
;				used then it is in units of the data range)
;
;	range=range	scalar, array[ndim] or array[2,ndim]
;				if scalar then the output array is scaled to the range
;					[0,range] in every dimension
;				if array[ndim] then the output array is scaled to the ranges
;					[0,range[i]], i=0,ndim-1, in each dimension
;				if array[2,ndim] then the output array is scaled to the ranges
;					[range[0,i],range[1,i]], i=0,ndim-1 in each dimension
;
;	/multid		(only in ndim > 1) reforms output array to [ndim,n] array.
;				The default is a two-dim array of size
;					[ndim,n[0]*n[1]*..*n[ndim-1]]
; OUTPUTS:
;	Result		array[ndim,n[0],n[1],..,n[ndim-1]]; type: long integer or float
;					coordinates across the area
;					a float array is returned if keywords /one or /range
;					are used.
; CALLS:
;	InitVar, BadValue, gridgen1d
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; RESTRICTIONS:
;	For time structures only 1-dim arrays can be generated
; PROCEDURE:
; >	The origin is subtracted after applying the scaling implied by
;	keywords /one and range=range
; >	Combination of 'replicate' function and matrix multiplication #
;
; >	Create 1-dim time grids:
;	gridgen(10,TimeUnit(/day))
;		10-element grid of time differences with 1-day time steps
;	gridgen(10,TimeUnit(/day),range='2008_001'
;	gridgen(10,range=['2008_001','2008_009'])
;		10-element grid with 1-day steps starting at 2008_001
; MODIFICATION HISTORY:
;	OCT-2000, Paul Hick (UCSD/CASS)
;		Result of merging indgen1d, indgen2d, and indgen3d
;	JAN-2002, Paul Hick (UCSD/CASS)
;		Change all replicate(1.0,n) statements to replicate(1L,n) in the
;		last case block.
;		Using 1.0 will always return a float array even when an integer
;		return would be acceptable.
;	NOV-2005, Paul Hick (UCSD/CASS)
;		Added 4D version
;	SEP-2006, Paul Hick (UCSD/CASS)
;		Substantial rewrite: generalized to work for any number
;		of dimensions (memory permitting)
;	OCT-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added argument 'unit' and allowed keywords 'origin' and
;		'range' to be standard structures. This allows setting up
;		1-dim arrays of times.
;-

InitVar, multid, /key

; This is a bit sick, but I have done worse

InitVar, merge, /key

IF merge THEN BEGIN

	CASE n_params()	OF
	1: r = nn
	2: BEGIN
		n1 = n_elements(nn)
		n2 = n_elements(yy)

		r = make_array(dim=[2,n1*n2], type=(IsType(nn) > IsType(yy)), /nozero)

 		r[0,*] = SuperArray(nn,n2,/trail)
		r[1,*] = SuperArray(yy,n1,/lead )

		IF multid THEN r = reform(r,2,n1,n2,/overwrite) ELSE r = reform(r,/overwrite)
	END
	3: BEGIN
		n1 = n_elements(nn)
		n2 = n_elements(yy)
		n3 = n_elements(zz)

		r = make_array(dim=[3,n1*n2*n3], type=(IsType(nn) > IsType(yy) > IsType(zz)), /nozero)

 		r[0,*] = SuperArray(nn,n2*n3,/trail)
		r[1,*] = SuperArray(SuperArray(yy,n1,/lead),n3,/trail)
		r[2,*] = SuperArray(zz,n1*n2,/lead)

		IF multid THEN r = reform(r,3,n1,n2,n3,/overwrite) ELSE r = reform(r,/overwrite)
	END
	ELSE: message, '/merge can only handle 1-,2- and 3-D grids'
	ENDCASE

	RETURN, r

ENDIF

IF IsType(range ,/string) THEN range  = TimeSet(range )
IF IsType(origin,/string) THEN origin = TimeSet(origin)
time_grid = IsType(unit,/defined) OR IsTime(range) OR IsTime(origin)

InitVar, nn, 0			; nn may be absent
n = long(nn)

ndim = n_elements(n)
multid AND= ndim GT 1

has_origin	= IsType(origin , /defined)
has_edge	= IsType(edge	, /defined)
has_open	= IsType(open	, /defined)
has_one 	= IsType(one	, /defined)

has_range = IsType(range, /defined)
IF has_range THEN BEGIN

	range_ = range
	IF time_grid AND n[0] EQ 0 THEN BEGIN
		InitVar, unit, TimeUnit(/day)
		n = round(TimeOp(/subtract,range_[1],range_[0],unit))+1
		range_ = range_[0]
	ENDIF

	CASE n_elements(range_) OF
	1   : BEGIN
		CASE time_grid OF
		0: range_ = [0,range_]
		1: BEGIN
			InitVar, unit, TimeUnit(/day)
			range_ = [range_,TimeOp(/add,range_,TimeSet(n-1,unit,/diff))]
		END
		ENDCASE
	END
	ndim: range_ = [intarr(1,ndim), reform(range_,1,ndim)]
	ELSE: range_ = reform(range_,2,ndim)
	ENDCASE

ENDIF

n_origin = n_elements(origin)  -1
n_edge	 = n_elements(edge	)  -1
n_open	 = n_elements(open	)  -1
n_one	 = n_elements(one	)  -1
n_range  = n_elements(range_)/2-1

; The original gridgen went upto ndim=4 using the following lines:
; r[0,*] = reform( reform( gridgen1d(n[0],origin=origin[0],edge=edge[0],open=open[0],one=one[0],range=range_[0:1],double=double) # replicate(1L,n[1]), n01 ) # replicate(1L,n[2]), n012 ) # replicate(1L,n[3])
; r[1,*] = reform( reform( replicate(1L,n[0]) # gridgen1d(n[1],origin=origin[1],edge=edge[1],open=open[1],one=one[1],range=range_[2:3],double=double), n01 ) # replicate(1L,n[2]), n012 ) # replicate(1L,n[3])
; r[2,*] = reform( reform( replicate(1L,n[0]) # replicate(1L,n[1]), n01 ) # gridgen1d(n[2],origin=origin[2],edge=edge[2],open=open[2],one=one[2],range=range_[4:5],double=double), n012 ) # replicate(1L,n[3])
; r[3,*] = reform( reform( replicate(1L,n[0]) # replicate(1L,n[1]), n01 ) # replicate(1L,n[2]), n012 ) # gridgen1d(n[3],origin=origin[3],edge=edge[3],open=open[3],one=one[3],range=range_[6:7],double=double)

; Calculate number of elements in grid
; Remember that gridgen1d adds an element if keyword edge is set.

np = lonarr(ndim)
IF has_edge THEN np[0] = n[0]+edge[0] ELSE np[0] = n[0]
FOR i=1,ndim-1 DO IF has_edge THEN np[i] = np[i-1]*(n[i]+edge[i < n_edge]) ELSE np[i] = np[i-1]*n[i]

FOR i=0,ndim-1 DO BEGIN 			; Loop over all ndim components

	FOR j=0,i-1 DO BEGIN
		grid = replicate(1L,n[j])
		IF j EQ 0 THEN tmp = grid ELSE tmp = reform( tmp # grid, np[j] )
	ENDFOR

	IF has_origin THEN origin_i	= origin[  i < n_origin]
	IF has_edge   THEN edge_i 	= edge  [  i < n_edge  ]
	IF has_open   THEN open_i 	= open	[  i < n_open  ]
	IF has_one	  THEN one_i	= one	[  i < n_one   ]
	IF has_range  THEN range_i	= range_[*,i < n_range ]

	grid = gridgen1d(n[i],unit,origin=origin_i,edge=edge_i,open=open_i,one=one_i,range=range_i,double=double)

	CASE i NE 0 OF
	0: BEGIN
		tmp = grid
		; The first call to gridgen1d sets the type of the output array.
		CASE time_grid OF
		0: r = make_array(type=size(grid,/type), dimension=[ndim,np[ndim-1]])
		1: r = replicate( !TheTime, [ndim,np[ndim-1]] )
		ENDCASE
	ENDCASE
	1: tmp = reform(tmp # grid, np[i])
	ENDCASE

	FOR j=i+1,ndim-1 DO BEGIN
		grid = replicate(1L,n[j])
		tmp = reform( tmp # grid, np[j] )
	ENDFOR

	r[i,*] = tmp

ENDFOR

IF multid THEN r = reform(r, [ndim,n], /overwrite) ELSE r = reform(r, /overwrite)

RETURN, r  &  END