;+
;
; NAME:
;
;       FT_SURFACE
;
;
; PURPOSE:
;
;		Computes the FFT of a real function and its right frequency axes scaling, labels and units.
;
; DESCRIPTION:
;
;		The FFT is always a bit difficult to graph: man's mind feels not easy to think in Fourier space,
;		so that each time you must think about what means the frequency axis of your FFT plot with
;		respect to your x units. More the IDL FFT vector contains both positive and negative frequencies
;		organized in a complicated format: first nu=0, then +nu_min to +nu_max then again -nu_max to -nu_min.
;		This routine made for you the correct surface plot of your FFT giving an easy to
;		understand graph of it.
;
;		FT_SURFACE computes the value and the right frequency axes scaling, labels and units of the
;		Fourier Transform of a real function, and draws a 2D shaded surface.
;		The user can choose to draw the Modulus, the Squared Modulus, the Real or Imaginary part or
;		the Phase of the FT, or just to make the computations without graphic outputs.
;		The input function z(x,y), i.e. the argument of the Fourier Tranform, can also be drawn
;		if x,y and z are supplied.
;
;
; CALLING SEQUENCE:
;
;		FT_Surface, FT, FT=FT, DATA_X=DATA_X, DATA_Y=DATA_Y, DATA_Z=DATA_Z, MODE=MODE, /NORMALIZE, $
;				XY_UNIT=XY_UNIT, NU_FACTOR=NU_FACTOR, NU_TITLE=NU_TITLE, $
;				NU_X_AXIS=NU_X_AXIS, NU_Y_AXIS=NU_Y_AXIS, OUT_ZLABEL=OUT_ZLABEL, OUTPUTZ=OUTPUTZ, $
;				/XY_ISOMETRIC, /SHOW_DATA, /LOG_DATA, DATA_TITLE=DATA_TITLE, /NO_PLOT
;
; INPUTS:
;
;       FT:		(obsolete) The Fourier Transform to be plotted: 2D array of real or complex data.
;				Must be computed before to call FT_SURFACE, e.g. by means of the FFT function.
;
;
; KEYWORDS:
;
;		FT:		If the parameter FT is not specified (considered obsolete) the FT is computed internally
;				and the FT keyword can be used to output it.
;				NOTE: the FT keyword always returns the FFT of the whole DATA_Z array, even if CROP_FACTOR is used.
;
;		DATA_X:	Vector of X coordinates of the original data z(X,Y) of which FT is the transform.
;				If not supplied, an integer vector from 0 to x size of z is used for x-axis.
;
;		DATA_X:	Same meaning of DATA_X but for y-axis.
;
;		DATA_Z:	2D array Z values of the original data z(X,Y) of which FT is the transform.
;				If not supplied, the keyword /SHOW_DATA cannot be set.
;
;		MODE:	String that specifies the type of plot as follows:
;				'SQUARED MODULUS': 	the default, plot the squared modulus of FT;
;				'MODULUS':			plot the modulus of FT;
;				'REAL':				plot the real part of FT;
;				'IMAGINARY':		plot the imaginary part of FT;
;				'PHASE':			plot the phase of FT;
;				'COMPLEX':			output the complex FT, NO_PLOT keyword is set;
;
;		CROP_FACTOR: tells the program that the output is desired only until a maximum frequency which is CROP_FACTOR
;				times the maximum input frequency.
;				This is useful when not aliased frequency-limited spectrum is required for a function which has non zero power
;				at frequencies higher than this limit: in this case the only way to compute a not-aliased result is to input
;				an highly sampled DATA_Z array and take only its spectrum within the frequency limit.
;				NOTE: if the CROP_FACTOR value does not led to an integer number of pixels, the limit frequency is not reached
;				and the largest integer less than the fractional value is used.

;		NORMALIZE: if set, the plotted modulus or squared modulus of FT are normalized to the maximum.
;				Does not work in 'real', 'imaginary' or 'phase' modes.
;
;		XY_UNIT: String: units of measure of the X and Y coordinates of the original data Z(X,Y) of which FT
;				is the transform. We suggest to always supply the units of measure of X and Y,
;				so that the routine	can compute the correct scaling of the frequency nu-axis
;				in units of cycles/XY_UNIT.
;				Does not work if DATA_X or DATA_Y are not supplied.
;
;		NU_FACTOR: Float: if you want to draw the FT versus a custom nu-axis with units different
;				than cycles/XY_UNIT, just specify here the multiplying factor:
;				your_nu_unit = NU_FACTOR * (cycles/XY_UNIT).
;				If this keyword is set, we suggest to also set NU_TITLE to the final frequency label and unit.
;
;		NU_TITLE: String: label of the frequency axis with units, useful when NU_FACTOR is different from 1.
;
;		NU_X_AXIS: Output array with the right nu-x_axis: the nu-x_axis is defined as follows:
;
;				Even X-size of FT:
;
;					-nu_max, -nu_(max-1), .... , -nu_min, 0, nu_min, ... , nu_(max-1)
;
;				Odd X-size of FT:
;
;					-nu_max, .... , -nu_min, 0, nu_min, ... , nu_max
;
;				NOTE: 	if the data has an even number of elements, then the frequencies goes from
;						zero to the Nyquist frequency nu_max = N/2 * nu_factor / data_range, else if it has
;						an odd number of points, then the Nyquist frequency is not evaluated, because
;						such data are completely determined from the frequencies from zero to
;						nu_max = (N/2-.5) * nu_factor / data_range.
;
;		NU_Y_AXIS: output array with the nu-y_axis: the nu-y_axis is defined the same as the nu-x_axis.
;
;		OUTPUTZ: output 2D-array with the required mode of power spectra.
;
;		OUT_ZLABEL: output the string of the z-axis label with units.
;
;		SHOW_DATA: if set, the original data Z(X,Y) is also drawn. DATA_Z must be supplied.
;				If DATA_Z is a COMPLEX array, then its MODULUS is displayed.
;				If SHOW_DATA is an integer greater than 31, it is assumed to be the window number of an
;				existing Draw_Widget.
;
;		LOG_DATA: if set and also SHOW_DATA is set, the data Z(X,Y) are plotted in a /zlog graph.
;
;		DATA_TITLE: if SHOW_DATA is set, title string for the data surface plot.
;
;		NO_PLOT: if set, then the FT surface is not plotted: use to only get the output data.
;				The original data surface Z(X,Y) can still be plotted if /SHOW_DATA is set.
;
;		_EXTRA:	Any keyword of the SHADE_SURFACE routine can be used in calling this function.
;				Caution in using the keywords xtitle, ytitle and ztitle which will overwrite
;				the automatic titles and the NU_TITLE settings.
;
;
; RESTRICTIONS:
;
;		This routine is only suitable for use with FFT of *real functions*, not for FFT of complex data.
;
;
; MODIFICATION HISTORY:
;
;       September 2004, G. Li Causi, Rome Astronomical Observatory
;						licausi@mporzio.astro.it
;						http://www.mporzio.astro.it/~licausi/
;		May 2005: G.Li Causi - SHOW_DATA can now receive a draw widget window number
;							 - DATA_TITLE added
;							 - NO_PLOT added
;		September 2005: G.Li Causi - COMPLEX mode added
;		October 2007: G Li Causi - Correct an error in frequency axis computation
;		May 2008: G. Li Causi - added CROP_FACTOR for frequency-limited output;
;								also work with complex DATA_Z input;
;								added FT keyword.
;-


PRO FT_Surface, FT0, FT=FT, DATA_X=DATA_X, DATA_Y=DATA_Y, DATA_Z=DATA_Z, CROP_FACTOR=CROP_FACTOR, $
			MODE=MODE, NORMALIZE=NORMALIZE, $
			XY_UNIT=XY_UNIT, NU_FACTOR=NU_FACTOR, NU_TITLE=NU_TITLE, DATA_TITLE=DATA_TITLE, $
			SHOW_DATA=SHOW_DATA, LOG_DATA=LOG_DATA, XY_ISOMETRIC=XY_ISOMETRIC, NO_PLOT=NO_PLOT, $
			NU_X_AXIS=NU_X_AXIS, NU_Y_AXIS=NU_Y_AXIS, OUTPUTZ=OUTPUTZ, OUT_ZLABEL=OUT_ZLABEL, $
            _EXTRA=EXTRA


;************
;Input Check:
;************
;IF n_elements(ft0) EQ 0 THEN Message, 'No FT data has been provided!'
IF n_elements(ft0) EQ 0 THEN ft0 = FFT(DATA_Z, -1, /DOUBLE)

ft = ft0

s = size(ft)
IF s[0] NE 2 THEN Message, 'FT array must be 2-dimensional!'

IF NOT KEYWORD_SET(data_x) THEN data_x = dindgen(s[1])
IF NOT KEYWORD_SET(data_y) THEN data_y = dindgen(s[2])

IF n_elements(data_x) NE s[1] THEN Message, 'X and Y vectors must be have same size of FT array!'
IF n_elements(data_y) NE s[2] THEN Message, 'X and Y vectors must be have same size of FT array!'

IF NOT KEYWORD_SET(NU_FACTOR) THEN NU_FACTOR = 1.

IF KEYWORD_SET(XY_UNIT) THEN data_unit = XY_UNIT ELSE data_unit = 'X'

IF NOT KEYWORD_SET(MODE) THEN MODE = 'SQUARED MODULUS'

IF NOT KEYWORD_SET(CROP_FACTOR) THEN CROP_FACTOR = 1
CROP_FACTOR = ABS(CROP_FACTOR)

MODE = strupcase(MODE)



;**********************
;Show data if required:
;**********************
win = !d.window
IF KEYWORD_SET(SHOW_DATA) AND NOT KEYWORD_SET(NO_PLOT) THEN BEGIN

	IF NOT KEYWORD_SET(data_z) THEN Message, 'No Z array have been supplied!'

	sz = size(data_z)
	IF sz[0] NE 2 THEN Message, 'Z array must be 2-dimensional!'
	IF sz[1] NE s[1] OR sz[2] NE s[2] THEN Message, 'Z array must the same size as FT!'

	IF SHOW_DATA GT 31 THEN BEGIN			;widget window
		wset, show_data
	ENDIF ELSE BEGIN
		IF win GT 0 THEN window, 0, title='Input Data' ELSE window, 1, title='Input Data'
	ENDELSE


	IF SIZE(data_z, /TNAME) EQ 'COMPLEX' THEN BEGIN
		data_z_plot = ABS(data_z)
		z_title = 'Modulus'
	ENDIF ELSE BEGIN
		data_z_plot = data_z
		z_title = 'Intensity'
	ENDELSE

	zmin = min(data_z_plot)
	IF KEYWORD_SET(LOG_DATA) THEN BEGIN
		zlog=1
		IF zmin LE 0 THEN zmin = min(data_z_plot[where(data_z_plot GT 0)])
	ENDIF

	IF KEYWORD_SET(XY_ISOMETRIC) THEN BEGIN
		x_range = [min(data_x) < min(data_y), max(data_x) > max(data_y)]
		y_range = x_range
	ENDIF ELSE BEGIN
		x_range = [min(data_x), max(data_x)]
		y_range = [min(data_y), max(data_y)]
	ENDELSE

	index = where(data_z_plot GE zmin, count)
	IF count EQ 1 THEN BEGIN
		zlog = 0
		zmin = min(data_z_plot)
	ENDIF

	;Plot Surface:
	SHADE_SURF, data_z_plot, data_x, data_y, zlog=zlog, xtitle=data_unit, ytitle=data_unit, zrange=[zmin, max(data_z_plot)], $
			/xstyle, /ystyle, /zstyle, ztitle=z_title, charsize=3, yrange=y_range, xrange=x_range, MIN_VALUE=zmin

	;Main title:
	IF KEYWORD_SET(DATA_TITLE) THEN	XYOUTS, 0.5, 0.95, DATA_TITLE, ALIGNMENT=0.5, /NORMAL

ENDIF



;*******************************************
;Compute input nu_x and nu_y frequency axes:
;*******************************************
n_samples_x = s[1]
n_samples_y = s[2]

;Center the origin:
ft = SHIFT(ft, n_samples_x/2, n_samples_y/2)

;Frequency ranges definition:
n_nu_points_x = n_samples_x				;numero elementi parte positiva e negativa frequenze FT
n_nu_points_y = n_samples_y				;numero elementi parte positiva e negativa frequenze FT

x_range = max(data_x) - min(data_x) 	;total x range of the data
y_range = max(data_y) - min(data_y) 	;total y range of the data

x_parity = FIX(n_nu_points_x/2.) NE (n_nu_points_x/2.)		;0=even, 1=odd
y_parity = FIX(n_nu_points_y/2.) NE (n_nu_points_y/2.)

;;;nu_x = (dindgen(n_nu_points_x)/(n_nu_points_x - x_parity) - 0.5) / (x_range / n_samples_x)	;x-axis in cycles/XY_UNITS	<--sbagliato, infatti per n_points dispari NON si raggiunge la freq di nyquist 0.5
;;;nu_y = (dindgen(n_nu_points_y)/(n_nu_points_y - y_parity) - 0.5) / (y_range / n_samples_y)	;y-axis in cycles/XY_UNITS

nu_x = [ dindgen(n_nu_points_x/2 + 1) , -reverse(dindgen(n_nu_points_x/2 - 1 + x_parity) + 1) ] / n_samples_x * ((n_samples_x-1) / x_range)	;nu_x-axis in cycles/XY_UNITS
nu_y = [ dindgen(n_nu_points_y/2 + 1) , -reverse(dindgen(n_nu_points_y/2 - 1 + y_parity) + 1) ] / n_samples_y * ((n_samples_y-1) / x_range)	;nu_y-axis in cycles/XY_UNITS

nu_x = SHIFT(nu_x, - (n_nu_points_x/2 + 1))
nu_y = SHIFT(nu_y, - (n_nu_points_y/2 + 1))

nu_x = nu_x * NU_FACTOR		;nu_x-axis in user units: NU_FACTOR * (cycles/XY_UNIT)
nu_y = nu_y * NU_FACTOR		;nu_y-axis in user_units: NU_FACTOR * (cycles/XY_UNIT)


;**********************
;Output frequency axes:
;**********************
out_n_samples_x = FIX(n_samples_x * CROP_FACTOR)
out_n_samples_y = FIX(n_samples_y * CROP_FACTOR)

NU_X_AXIS = nu_x[n_samples_x/2. - out_n_samples_x/2.: n_samples_x/2. + out_n_samples_x/2. - 1]
NU_Y_AXIS = nu_y[n_samples_y/2. - out_n_samples_y/2.: n_samples_y/2. + out_n_samples_y/2. - 1]

nu_x_range = x_range / NU_FACTOR * CROP_FACTOR		;nu x_range in user_units
nu_y_range = y_range / NU_FACTOR * CROP_FACTOR		;nu y_range in user_units

;plot_x_range = [-out_n_samples_x/nu_x_range/2, out_n_samples_x/nu_x_range/2] * CROP_FACTOR
;plot_y_range = [-out_n_samples_y/nu_y_range/2, out_n_samples_y/nu_y_range/2] * CROP_FACTOR

plot_x_range = minmax(nu_x_axis)
plot_y_range = minmax(nu_y_axis)


IF KEYWORD_SET(XY_ISOMETRIC) THEN BEGIN
	plot_x_range = [min(plot_x_range[0]) < min(plot_y_range[0]), max(plot_x_range[1]) > max(plot_y_range[1])]
	plot_y_range = plot_x_range
ENDIF


;Axis units and Titles:
IF NOT KEYWORD_SET(XY_UNIT) THEN XY_UNIT = 'units of X'

IF NU_FACTOR EQ 1 THEN nu_fac_string = '' ELSE nu_fac_string = 'units of ' + strtrim(string(NU_FACTOR),2)

IF KEYWORD_SET(NU_TITLE) THEN xtitle = NU_TITLE $
	ELSE xtitle='Frequency (' + nu_fac_string + ' cycles/' + XY_UNIT + ')'

ytitle = xtitle
ztitle = ''



;***************
;MODE Selection:
;***************
CASE MODE OF

	'REAL': BEGIN
		OUTPUTZ = FLOAT(ft[n_samples_x/2. - out_n_samples_x/2.: n_samples_x/2. + out_n_samples_x/2. - 1, n_samples_y/2. - out_n_samples_y/2.: n_samples_y/2. + out_n_samples_y/2. - 1])
		zrange = [min(OUTPUTZ), max(OUTPUTZ)]
		ztitle = ztitle + ' Real Part'
	END

	'IMAGINARY': BEGIN
		OUTPUTZ = IMAGINARY(ft[n_samples_x/2. - out_n_samples_x/2.: n_samples_x/2. + out_n_samples_x/2. - 1, n_samples_y/2. - out_n_samples_y/2.: n_samples_y/2. + out_n_samples_y/2. - 1])
		zrange = [min(OUTPUTZ), max(OUTPUTZ)]
		ztitle = ztitle + ' Imaginary Part'
	END

	'MODULUS': BEGIN
		OUTPUTZ = ABS(ft[n_samples_x/2. - out_n_samples_x/2.: n_samples_x/2. + out_n_samples_x/2. - 1, n_samples_y/2. - out_n_samples_y/2.: n_samples_y/2. + out_n_samples_y/2. - 1])
		zrange = [0, max(OUTPUTZ)]
		ztitle = ztitle + ' Modulus'
	END

	'SQUARED MODULUS': BEGIN
		OUTPUTZ = (ABS(ft[n_samples_x/2. - out_n_samples_x/2.: n_samples_x/2. + out_n_samples_x/2. - 1, n_samples_y/2. - out_n_samples_y/2.: n_samples_y/2. + out_n_samples_y/2. - 1]))^2
		zrange = [0, max(OUTPUTZ)]
		ztitle = ztitle + ' Squared Modulus'
	END

	'PHASE': BEGIN
		OUTPUTZ = ATAN(ft[n_samples_x/2. - out_n_samples_x/2.: n_samples_x/2. + out_n_samples_x/2. - 1, n_samples_y/2. - out_n_samples_y/2.: n_samples_y/2. + out_n_samples_y/2. - 1], /PHASE)
		zrange=[-!PI, !PI]
		ztitle = 'Phase (radians)'
	END

	'COMPLEX': BEGIN
		OUTPUTZ = ft[n_samples_x/2. - out_n_samples_x/2.: n_samples_x/2. + out_n_samples_x/2. - 1, n_samples_y/2. - out_n_samples_y/2.: n_samples_y/2. + out_n_samples_y/2. - 1]
		NO_PLOT = 1			;set no plotting for complex mode
	END

ENDCASE


;**************
;Normalization:
;**************
IF KEYWORD_SET(NORMALIZE) AND (MODE EQ 'MODULUS' OR MODE EQ 'SQUARED MODULUS') THEN BEGIN
	OUTPUTZ = OUTPUTZ / MAX(OUTPUTZ)
	zrange = [0, 1]
	ztitle = 'Normalized ' + ztitle
ENDIF

OUT_ZLABEL = ztitle


;*********
;Plotting:
;*********
IF NOT KEYWORD_SET(NO_PLOT)THEN BEGIN

	IF win GT 0 THEN BEGIN
		wset, win
		wshow, win
	ENDIF ELSE window, 0, title='Fourier Transform'

	;Title:
	IF n_elements(EXTRA) NE 0 THEN BEGIN
		IF KEYWORD_SET(EXTRA.TITLE) THEN BEGIN
			TITLE = EXTRA.TITLE
			EXTRA.TITLE = ''
		ENDIF ELSE TITLE = ''
	ENDIF

	;Plot Surface:
	SHADE_SURF, OUTPUTZ, NU_X_AXIS, NU_Y_AXIS, $
		zrange=zrange, xrange=plot_x_range, yrange=plot_x_range, $
		/xstyle, /ystyle, /zstyle, $
		xtitle=xtitle, ytitle=ytitle, ztitle=ztitle, TITLE='', $
		_EXTRA=EXTRA

	;Main title:
	IF n_elements(TITLE) NE 0 THEN XYOUTS, 0.5, 0.95, TITLE, ALIGNMENT=0.5, /NORMAL

ENDIF


END