first commit slope_computation_package folder

;+ Author: Laura Schreiber - August 2022

;

;ROUTINE NAME:

;asong_leasy_tilt

;

;PURPOSE:

;Computes the tilt signal x and y from the signal in each pupil using the formula:

;tilt x = ((A+B)-(C+D))/ A+B+C+D;      tilt y = ((A+C)-(B+D))/ A+B+C+D;      

;

;CALLING SEQUENCE:

;asong_leasy_tilt, frame, qmask, tilt

;

;INPUT:

;frame:  2D (3D if SEQ keyword is set) array containing the image of the 4 pupil 

;qmask:  3D array with four planes, each representing a binary mask 

;        associated to one quadrant of the sensor

;        

;KEYWORDS:

;SEQ: in case of sequential pupil images (no MLA to produce 4 pupil images)  

;GAIN : gain to be computed or calibrated extenally        

;

;OUTPUT: 

;tilt :  vector of 2 elements containin the tilt x and y signal 

;-

; NOTE: the tilt is in signal unless the gain is passed to the procedure

PRO asong_leasy_tilt, frame, qmask, $ ; inputs

                      SEQ = SEQ,  GAIN = GAIN,  $ ; keywords

                      tilt            ; output

if keyword_set(SEQ) then begin

  A = total(qmask[*,*,0]*frame[*,*,0])

  B = total(qmask[*,*,1]*frame[*,*,1])

  C = total(qmask[*,*,2]*frame[*,*,2])

  D = total(qmask[*,*,3]*frame[*,*,3])  

endif else begin  

  A = total(qmask[*,*,0]*frame)

  B = total(qmask[*,*,1]*frame)

  C = total(qmask[*,*,2]*frame)

  D = total(qmask[*,*,3]*frame)

endelse

tot = A+B+C+D

;tilt=[((A+B)-(C+D))/tot,((A+C)-(B+D))/tot]

tilt=[((B+D)-(A+C))/tot,((C+D)-(A+B))/tot]

;tilt=[((B+C)-(A+D))/tot,((C+D)-(A+B))/tot]

RETURN

END

 No newline at end of file

;+

; NAME:

;  pwd

; PURPOSE:

;  Print current working directory.

; DESCRIPTION:

;

; CATEGORY:

;  Miscellaneous

; CALLING SEQUENCE:

;  pwd

; INPUTS:

;  None.

; OPTIONAL INPUT PARAMETERS:

;

; KEYWORD INPUT PARAMETERS:

;

; OUTPUTS:

;

; KEYWORD OUTPUT PARAMETERS:

;

; COMMON BLOCKS:

;

; SIDE EFFECTS:

;

; RESTRICTIONS:

;

; PROCEDURE:

;

; MODIFICATION HISTORY:

;  1999 May 13, Written by Marc W. Buie, Lowell Observatory

;

;-

pro pwd

   cd,current=current

   print,current

end

;+ Author: Laura Schreiber - August 2022

;

; ROUTINE NAME:

;SEQ_asong_slope

;

;PURPOSE:

;Given a frame and the binary masks defining the four quadrants of the ASONG WFS, compute the wavefront derivatives in X and Y.

;Modification of the asong_slope procedure, when the images of the 4 pupils are taken in sequence. 

;

;CALLING SEQUENCE:

;asong_slope, frame, qmask, pup, centres, sx, sy

;

;INPUT:

;frame:  cube of measured frames with subtracted background and divided by the flat (if required)

;qmask:  3D array with four planes, each representing a binary mask 

;        associated to one quadrant of the sensor

;pup:    footprint of sx and sy

;centres: 4 X 2-elements vector of coordinates of pupils centroids

;

;KEYWORDS:

;gain:        when this keyword is set, the gain is not computed and the value set in the keyword is used instead.

;             It can come from calibration procedures. DEFAULT VALUE = 0

;pyr:         Set this keyword in case a Pyramid WFS is used in stead of ASONG. DEFAULT VALUE = 0    

;dxf, dyf:    GTF x and y period in mm. If only one of the two keywords is set, the period in X and Y is considered to be the same. 

;             DEFAULT VALUE = 9.05 mm  

;foc:         distance from the aperture and the fourier plane

;

;OUTPUT:

;sx, sy: array of signals in X and Y

;-

PRO SEQ_asong_slope, frame, qmask, pup, centres, $                            ;inputs

                  GAIN = GAIN, dxf = dxf, dyf = dyf, foc = foc, $             ;keywords

                  sx, sy                                                      ;outputs

;on_error, 2

if not keyword_set(PYR) then begin

   if not keyword_set(foc) then foc = 500.  ; focal lengh in mm

   if not keyword_set(dxf) then dxf = 9.05    ; period in mm

   dyf = dxf

   gx = dxf / foc / 2. / !dpi / sqrt(2.)

   gy = dyf / foc / 2. / !dpi / sqrt(2.)

endif else begin 

   if not keyword_set(GAIN) then begin

      gx = 1.0 & gy = 1.0  

   endif else begin

      gx = gain & gy = gain

   endelse

endelse

dimpup = (size(pup,/dim))[0]

c=centres

c = transpose(c)

; Extract four quadrants from the frame

nquad = 4

quad = fltarr(dimpup, dimpup, nquad)

lx = round(c[*,0]) - dimpup/2

ux = round(c[*,0]) + dimpup/2 

ly = round(c[*,1]) - dimpup/2

uy = round(c[*,1]) + dimpup/2

for q = 0, nquad - 1 do begin

   quad[*,*,q] = (frame[*,*,q] * qmask[*,*,q])[lx:ux-1,ly:uy-1]*pup

endfor

; Compute signal and derivative

sx = fltarr(dimpup, dimpup)

sy = fltarr(dimpup, dimpup)

norm = total(quad, 3)

quadmask = (pup and (norm ne 0)) and 1B

w = where(quadmask)

if keyword_set(PYR) then begin

  sx[w] = ((quad[*,*,3] + quad[*,*,1]) - (quad[*,*,2] + quad[*,*,0]))[w]

  sy[w] = ((quad[*,*,2] + quad[*,*,3]) - (quad[*,*,0] + quad[*,*,1]))[w]

  sx[w] = gx * sx[w] / norm[w]

  sy[w] = gy * sy[w] / norm[w]  

endif else begin

  sx[w] = ((quad[*,*,3] + quad[*,*,1]) - (quad[*,*,2] + quad[*,*,0]))[w]/visx

  sy[w] = ((quad[*,*,0] + quad[*,*,1]) - (quad[*,*,2] + quad[*,*,3]))[w]/visy

  sx[w] = sx[w] / norm[w]

  sy[w] = sy[w] / norm[w]  

  sx[w] = asin(sx[w]) * gx

  sy[w] = asin(sy[w]) * gy  

endelse

return

end