first commit wfe_reconstruction_packega folder

; $Id: array_overlap.pro, v 1.0 Aug 1999 e.d. $

;

;+

; NAME:

;	ARRAY_OVERLAP

;

; PURPOSE:

;	Find bounds of overlap region of two 2D arrays, by ideally

;	matching the coordinates of a reference point.

;

; CATEGORY:

;	Array manipulation.

;

; CALLING SEQUENCE:

;	ARRAY_OVERLAP, Size1, Size2, R1, R2, $

;   			   Lx1, Ux1, Ly1, Uy1, Lx2, Ux2, Ly2, Uy2

;

; INPUTS:

;	Size1:	2-components vector, size of array 1, as returned by

;		size52(array1, /DIM)

;

;	Size2:	2-components vector, size of array 2

;

;	R1:	2-components vector, coordinates of reference pixel in array 1

;

;	R2:	2-components vector, coordinates of reference pixel in array 2

;

; OUTPUTS:

;	Lx1, Ux1, Ly1, Uy1:	Lower and Upper X- and Y- bounds of intersection

;		in array 1

;

;	Lx2, Ux2, Ly2, Uy2:	Lower and Upper X- and Y- bounds of intersection

;		in array 2

;

; MODIFICATION HISTORY:

; 	Written by:	Emiliano Diolaiti, August 1999.

;-

; INT_BOUNDS: auxiliary procedure.

PRO int_bounds, p1, p2, size1, size2, l, u

	l = (p1 - p2) > 0 < (size1 - 1)

	u = (p1 + size2 - p2 - 1) > 0 < (size1 - 1)

	return

end

PRO array_overlap, size1, size2, r1, r2, $

   				   lx1, ux1, ly1, uy1, lx2, ux2, ly2, uy2

	on_error, 2

	int_bounds, r1[0], r2[0], size1[0], size2[0], lx1, ux1

	int_bounds, r1[1], r2[1], size1[1], size2[1], ly1, uy1

	int_bounds, r2[0], r1[0], size2[0], size1[0], lx2, ux2

	int_bounds, r2[1], r1[1], size2[1], size1[1], ly2, uy2

	return

end

;---------------------------------------------------------------------------------

;

;-Cette procedure reconstruit la WFE sur la pupille principale pas LS Fitting

;

;   Entrees : px          tableau des pentes en X estimates sur la pupille principale

;             py          tableau des pentes en Y estimates sur la pupille principale

;             masque      tableau des masques dans la pupille principale

;             recmat      tableau des recontruction (calcule par recmat.pro)

;             zer         zernike cube matrix

;

;   Sorties : wfe         tableau complexe de la WFE reconstruite sur la pupille principale

;

PRO calc_wfe_ls, px, py, masque, recmat, zer, wfe, c, notilt = notilt

  s = size(px,/dim)

  nzer = (size(recmat,/dim))[1]

  wfe = dblarr(s[0],s[1])

  c = transpose(recmat ## sort_subap(masque, px, py))

if not keyword_set(notilt) then $

  for i = 0,nzer-1 do wfe = wfe + c[i] * zer[*,*,i] else $

  for i = 2,nzer-1 do wfe = wfe + c[i] * zer[*,*,i]

return

end

 No newline at end of file

; $Id: coordinate.pro v 1.0 $

;

;+

; NAME:

; COORDINATE

;

; PURPOSE:

; This function defines a vector of coordinates symmetric with respect to an origin.

;

; CATEGORY:

; Array operations.

;

; CALLING SEQUENCE:

; Result = FUNCTION_NAME(N)

;

; INPUTS:

; N:  Integer number of elements. If N is not integer, it is truncated to 

;    the lower integer value.

;

; OUTPUTS:

; This function returns a float array of coordinates defined as follows:

; -(N-1)/2, ..., -0.5, +0.5, ..., +(N-1)/2, if N is even

; -N/2, ..., -1.0, 0.0, +1.0, ..., +N/2, if N is odd.

; 

; EXAMPLE:

; IDL> print, coordinate(4)

;     -1.50000    -0.500000     0.500000      1.50000

; IDL> print, coordinate(3)

;     -1.00000     0.000000      1.00000

;

; MODIFICATION HISTORY:

;   Written by: E.D. (INAF-OABo), April 2011

;-

function coordinate, n

k = fix(n)

x = findgen(k) - k/2

if k mod 2 eq 0 then x = x + 0.5

return, x

end

 No newline at end of file

; $Id: extend_array.pro, v 1.1 Jan 2000 e.d. $

;

;+

; NAME:

;	EXTEND_ARRAY

;

; PURPOSE:

;	Pad 2D array with 0s.

;

; CATEGORY:

;	Array manipulation.

;

; CALLING SEQUENCE:

;	Result = EXTEND_ARRAY(Array, S0, S1)

;

; INPUTS:

;	Array:	Array to be 0-padded

;

;	S0:	Size of extended array

;

; OPTIONAL INPUTS:

;	S1:	Second (y-) size of extended array. The default is S1 = S0

;

; KEYWORD PARAMETERS:

;	NO_OFF: 	Set this keyword to a nonzero value to indicate that the

;		input array must not be centered in the final frame. In this case

;		the element [0,0] of the output array coincides with the element

;		[0,0] of the input

;

;	POW2:	Set this keyword to a nonzero value to indicate that the array

;		size must be extended to the nearest power of 2. In this case, the

;		input values of S0 and S1 are overriden

;

; OUTPUTS:

;	Result:	0-padded array. If the keyword NO_OFF is not set, the original

;		array is centered in the final frame.

;		The input array is returned if the requested output size is smaller

;		than the input size

;

; OPTIONAL OUTPUTS:

;	OFFSET:	Use this output keyword to retrieve a 2-components vector with

;		the x- and y- offsets of the [0, 0] element of the input array

;		in the final frame

;

; SIDE EFFECTS:

;	If POW2 is set, the input values of S0 and S1 are overridden. If S0 and

;	S1 are named variables, their values are overwritten.

;

; RESTRICTIONS:

;	Apply only to 2D arrays.

;

; MODIFICATION HISTORY:

; 	Written by:	Emiliano Diolaiti, August 1999.

;	Updates:

;	1) Fixed bug on keyword POW2 (Emiliano Diolaiti, January 2000).

;	2) Fixed other bug on keyword POW2 (E. D., May 2014).

;-

FUNCTION extend_array, array, s0, s1, NO_OFF = no_off, POW2 = pow2, OFFSET = o

	on_error, 2

	if  size52(array, /N_DIM) ne 2  then  return, array

	if  n_params() eq 1 and not keyword_set(pow2)  then  return, array $

	else  if  n_params() eq 2  then  s1 = s0

	s = size52(array, /DIM)

	if  keyword_set(pow2)  then begin

	   s0 = s[0]  &  s1 = s[1]

	   l = log2(s0)  &  if  2^l ne s0  then  s0 = 2^(l + 1)

	   l = log2(s1)  &  if  2^l ne s1  then  s1 = 2^(l + 1)

	endif

	if  s0 eq s[0] and s1 eq s[1] or s0 lt s[0] or s1 lt s[1]  then begin

	   array1 = array  &  o = [0, 0]

	endif else begin

	   if  keyword_set(no_off)  then  o = [0, 0] $

	   else begin

	      o = [s0, s1] - s  &  o = (o + o mod 2) / 2

	   endelse

	   array1 = make_array(s0, s1, TYPE = size52(array, /TYPE))

	   array1[o[0],o[1]] = array

	endelse

	return, array1

end

; Copyright (c)  Harris Geospatial Solutions, Inc. All

;       rights reserved. Unauthorized reproduction is prohibited.

;+

; NAME:

;       FACTORIAL

;

; PURPOSE:

;       This function computes the factorial N! as the double-precision

;       product, (N) * (N-1) * (N-2) * ...... * 3 * 2 * 1.

;

; CATEGORY:

;       Special Functions.

;

; CALLING SEQUENCE:

;       Result = Factorial(n)

;

; INPUTS:

;       N:    A non-negative scalar or array of values.

;

; KEYWORD PARAMETERS:

;       STIRLING:    If set to a non-zero value, Stirling's asymptotic

;                    formula is used to approximate N!.

;

;       UL64: Set this keyword to return values as unsigned 64-bit integers.

;

; EXAMPLE:

;       Compute 20! with and without Stirling's asymptotic formula.

;         result_1 = factorial(20, /stirling)

;         result_2 = factorial(20)

;

;       Result_1 and result_2 should be 2.4227869e+18 and 2.4329020e+18

;       respectively.

;

; REFERENCE:

;       ADVANCED ENGINEERING MATHEMATICS (seventh edition)

;       Erwin Kreyszig

;       ISBN 0-471-55380-8

;

; MODIFICATION HISTORY:

;       Written by:  GGS, RSI, November 1994

;       CT, RSI, June 2000: Rewrote to handle array input; uses GAMMA for

;              non-integers; added UL64 keyword.

;       CT, RSI, October 2000: Separate calculation for scalar input.

;              Use common variable for look-up table.

;-

function factorial, x, stirling = stirling, UL64=ul64

  COMPILE_OPT idl2

  on_error, 2

  COMMON commonFactorial, factorialBuiltIn

	IF KEYWORD_SET(stirling) THEN BEGIN

		if MIN(x) lt 0 then $

			message, 'Values for N must be non-negative.'

		;Approximate N! using Stirling's formula.

		RETURN, SQRT(2.0d * !DPI * x) * (x / EXP(1.0d))^(x+0.0d)

	ENDIF

	IF (N_ELEMENTS(factorialBuiltIn) LT 1) THEN BEGIN

		factorialBuiltIn = $

			[1ull,1,2,6,24,120,720,5040,40320,362880,3628800, $

			39916800,479001600,6227020800,87178291200,1307674368000, $

			20922789888000,355687428096000,6402373705728000, $

			121645100408832000,2432902008176640000ull]

	ENDIF

	isInt = (x EQ FIX(x)) AND (x LE 20)

	n = N_ELEMENTS(x)

; For speed purposes, split the calculation into 1-element or array input

	IF (n EQ 1) THEN BEGIN   ; scalar or 1-element vector

		if x[0] lt 0 then $

			message, 'Values for N must be non-negative.'

		IF isInt[0] THEN BEGIN

		; do the integer factorials

			fact = factorialBuiltIn[x]

			IF NOT KEYWORD_SET(ul64) THEN fact = DOUBLE(fact)

		ENDIF ELSE BEGIN

		; do floating-point factorials using GAMMA(n+1)

		; the x*0 ensures that a 1-element vector remains a vector

			fact = x*0 + GAMMA(x+1d)

			IF KEYWORD_SET(ul64) THEN fact = ULONG64(fact)

		ENDELSE

	ENDIF ELSE BEGIN   ; array

		if MIN(x) lt 0 then $

			message, 'Values for N must be non-negative.'

		fact = KEYWORD_SET(ul64) ? ULON64ARR(n) : DBLARR(n)

		; first do all the integer factorials

		whereInt = WHERE(isInt,nInt, $

			COMPLEMENT=whereNotInt, NCOMPLEMENT=nNotInt)

		IF (nInt GT 0) THEN BEGIN

			; most common cases...

			fact[whereInt] = factorialBuiltIn[x[whereInt]]

		ENDIF

		; now do all the floating-point factorials using GAMMA(n+1)

		IF (nNotInt GT 0) THEN $

			fact[whereNotInt] = GAMMA(x[whereNotInt]+1d)

	ENDELSE

	RETURN, fact

end