; Reconstructor process that checks to see if a new reconstructor needs to be
; triggered and if so, triggers it.
;
; Set up as widget even though there is no GUI for consistency with other 
; programs and so one can be built later if needed
;
; Marcos van Dam, December 2005
; Modification history:
;   10Dec2007: If in stationary mode, look at aoprvert instead of obrt to trigger reconstructor
;
; Loop every second
; Trigger if
;   The recgo keyword has been set 
;   OR The angle has changed by more than 1 deg
;   AND The tip-tilt loop is closed
;   AND The median intensity is greater than 20 counts
;
; Subroutines
;   LGS_CENTROID_GAIN: calculate the centroid gain in each subaperture for the LGS
;   ADJUST_PARM_B: set the Bayesian regularization parameter using the median intensity
;   FINDHOLES: determine whether the reconstructor has holes and hence is bad
;   RECONSTRUCT_KL: create a Karhunen-Loeve modal reconstructor
;   RECONSTRUCT_BAYESIAN: find the Bayesian reconstructor
;   RECONSTRUCT: do the preliminaries for any reconstructor
;   TRIGGER_RECON: determine if a new reconstructor needs to be triggered
;   RECONPROC_EVENT: the event handler for reconproc
;   RECONPROC: initialize the reconstructor
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
;;
@kidl
PRO LGS_CENTROID_GAIN  
; calculates the centroid gain for each subaperture when using the LGS
; should also work for LGS AO on Keck I when it is implemented  
  COMPILE_OPT idl2
                                ; finds the centroid gain for every subaperture
  
  laser_spot=SHOW('ao.recmnlgs',error=error0,status=status0,/notrace)
  laser_spot_elong=SHOW('ao.recmxlgs',error=error1,status=status1,/notrace)
  sod_height=SHOW('ao.aofccsal',error=error2,status=status2,/notrace)
  prvert=SHOW('ao.aoprvert',error=error3,status=status3,/notrace) 
  elev=SHOW('dcs.el',error=error4,status=status4,/notrace)
  telescope=SHOW('ao.aotel',error=error5,status=status5,/notrace)
  
  IF status0+status1+status2+status3+status4+status5 NE 0 THEN BEGIN
     MESSAGE,/INFO,'Error reading keywords'
     PRINT,[error0,error1,error2,error3,error4,error5]
     RETURN
  ENDIF
  
  height=sod_height-4200.       ; sodium height from sea level minus the elevation of the telescope
  dl=0.562                       ; size of the lenslet

; need to reverse engineer to find the thickness
; from the smallest and elongated spots at zenith
; use approximate expressions
;elongx=dxv*thickness/height^2.*206265.*sin(elev)
  
  CASE telescope OF
     'Keck I': BEGIN
        baseline=5.
        pixel_size=1.5          ; ??? use true pixel size
        laser_loc_x=9.5
        laser_loc_y=9.5
     END
     
     'Keck II': BEGIN
        baseline=10.
        pixel_size=1.5          ; pixel size in arcsec
        
        laser_angle=(116.6-prvert) ; 0 is the top of the pupil between subapertures 300 and 301
; the angle rotates anti clockwise with increasing prvert angle
        laser_loc_x=sin(laser_angle*!Pi/180.)*10.5+9.5
        laser_loc_y=cos(laser_angle*!Pi/180.)*10.5+9.5
     END
  ENDCASE     
   
  charge_diffusion=0.5          ; FWHM in pixels
  charge_diffusion=charge_diffusion*pixel_size
  
  thickness=SQRT(laser_spot_elong^2.-laser_spot^2.)*height^2./(206265.*baseline)

; Need to calculate the elongation as fn of elevation
; The elongation goes as sin(elev)

; determine the location of the laser in subaperture units
; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
; define zero laser angle to be the top of the pupil


; find the distance to each subaperture
; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  dx=FLTARR(20,20)
  dy=dx

  for c1=0,19 do begin
     for c2=0,19 do begin
        dx[c1,c2]=ABS(c1-laser_loc_x)*dl
        dy[c1,c2]=ABS(c2-laser_loc_y)*dl
     endfor
  endfor

; load the subaperture map
  sub_map=INTARR(20,20)
  OPENR,lun,'/kroot/rel/ao/qfix/data/IdlParms/sub_ap_map.txt',/get_lun
  READF,lun,sub_map
  FREE_LUN,lun

  where_sub_map=WHERE(sub_map)

  dxv=dx[where_sub_map]
  dyv=dy[where_sub_map]

; calculate the elongation in arcsec
; ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

; approximate expressions, no longer used
;elongx=dxv*thickness/height^2.*206265.*sin(elev)
;elongy=dyv*thickness/height^2.*206265.*sin(elev)

; exact expressions
  b1=(height-thickness/2)/TAN(elev) ; baseline 1
  b2=(height+thickness/2)/TAN(elev) ; baseline 2

  elongx=(ATAN((height+thickness/2.)/(dxv+b2))-ATAN((height-thickness/2)/(dxv+b1)))*206265.
  elongy=(ATAN((height+thickness/2.)/(dyv+b2))-ATAN((height-thickness/2)/(dyv+b1)))*206265.

  fwhmvec=FLTARR(608)
  fwhmvec[INDGEN(304)*2]=SQRT(elongx^2.+laser_spot^2.+charge_diffusion^2)
  fwhmvec[INDGEN(304)*2+1]=SQRT(elongy^2.+laser_spot^2.+charge_diffusion^2)
    
; convert FWHM to centroid gain
  arcsec_per_cent=(SQRT(2.*!pi)*fwhmvec)/(2.*2.355)
  status=MODIFY('ao.wscngn',DOUBLE(arcsec_per_cent),error=error,/notrace)
  IF status NE 0 THEN BEGIN
     MESSAGE,/INFO,'Error writing ao.wscngn keyword'
     MESSAGE,/INFO,error
  ENDIF  
END

FUNCTION ADJUST_PARM_B, medianIntensity
  IF medianIntensity GT 500. THEN alpha=2.
  IF medianIntensity LE 500. AND medianIntensity GT 100. THEN alpha=5.
  IF medianIntensity LE 100. AND medianIntensity GT 70. THEN alpha=10.
  IF medianIntensity LE 70. THEN alpha=20.
  RETURN, alpha
END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION FINDHOLES, act
  ;;;
  ;;;Function to determine whether there is a hole in the reconstructor.
  ;;;A hole is defined as non-active subaperture(s) surrounded 
  ;;;by active subapertures.
  ;;;
  ;;;inputs:
  ;;;act- the 304 array of the active subapertures (1=active, 0 inactive)
  ;;;
  ;;;The algorithm searches every non-active subaperture to see what
  ;;;non-active subapertures it is connected to. If a non-active subaperture
  ;;;is not connected to any other non-active subaperture
  ;;;then it is a 1 element hole. Then the algorithm loops through the arrays 
  ;;;of what each zero is connected to, adding each connected zero together. If
  ;;;there are zeros not connected to the outside of the pupil or the central
  ;;;obscuration, then we have a hole.
  ;;;
  ;;;returns:
  ;;;0 if there is NO hole
  ;;;1 if there is a hole
  ;;;
  ;;;Richard Clare
  ;;;Created: 30 November 2006
  ;;;Last modified: 4 December 2006
  COMPILE_OPT idl2
  
  nr=20 ;;;the number of rows = number of columns
  
  ;;;load in the subaperture map
  openr,lun,'/kroot/rel/ao/qfix/data/IdlParms/sub_ap_map.txt',/get_lun
  sub_map=intarr(nr,nr)
  readf,lun,sub_map
  free_lun,lun   
  wheresubmap20=where(sub_map)
  
  ;;;add an extra row/col of zeros around the outside
  sub_map2=intarr(nr+2,nr+2)
  sub_map2[1:nr,1:nr]=sub_map
  wheresubmap=where(sub_map2)
  
  temp=fltarr(nr,nr)
  temp[wheresubmap20]=act
  recon_map=fltarr(nr+2,nr+2)  
  recon_map[1:nr,1:nr]=temp
  
  ;;;find the zeros within the pupil 
  zeroindex=where(recon_map eq 0)  
  zeros=fltarr(nr+2,nr+2)
  zeros[zeroindex]=1
  
  ;;;this is the array of what zeros each zero is connected to
  conindex=fltarr(n_elements(zeroindex),5)-1 ;;;5 as 4-connected plus self
  
  nr=nr+2 ;;;to account for the extra rows/cols I added around the outside
  
  hole=where(sub_map2 eq 0) ;;;starting point is ones that are always zero
  
  ;;;loop over each zero to see what it is connected to
  for i=0,n_elements(zeroindex)-1 do begin
     
     if total(zeroindex[i] eq hole) ne 1 then begin 
     ;;;only find connected ones to zeros inside the pupil   
        
     ;;;the array of connected indices
        conindex[i]=zeroindex[i] ;;;first is itself
        
     ;;;check each of its 4 neighbors provided not on the edge
        
        col=zeroindex[i] mod nr
        row=zeroindex[i]/nr
        
        j=1 ;;;j keeps track of how many are connected to that one
        
     ;;;check the one below
        if row ne 0 then begin
           if zeros[zeroindex[i]-nr] eq 1 then begin 
              conindex[i,j]=zeroindex[i]-nr
              j=j+1
           endif
        endif
        
     ;;;check the one to the left
        if col ne 0 then begin 
           if zeros[zeroindex[i]-1] eq 1 then begin 
              conindex[i,j]=zeroindex[i]-1
              j=j+1
           endif
        endif
        
     ;;;check the one to the right
        if col ne (nr-1) then begin
           if zeros[zeroindex[i]+1] eq 1 then begin 
              conindex[i,j]=zeroindex[i]+1
              j=j+1
           endif
        endif
        
     ;;;check the one above
        if row ne (nr-1) then begin
           if zeros[zeroindex[i]+nr] eq 1 then begin 
              conindex[i,j]=zeroindex[i]+nr
              j=j+1
           endif
        endif
        
     ;;;stop if there is a zero not connected to anything
        if j eq 1 then return, 1
     endif
     
  endfor
  
  ;;;need to loop ? times to ensure get them all
  ntimes=7
  
;;;now need to loop through the conindex to join up the connected ones  
  
  for n=0,ntimes-1 do begin      
     for i=1,n_elements(zeroindex)-1 do begin ;;;loop over each zero
        if total(zeroindex[i] eq hole) ne 1 then begin 
     ;;;only check connected ones to zeros inside the pupil              
           for j=0, total(conindex[i,*] gt 0)-1 do begin 
;;;loop over each one the zero is connected to              
              for k=0, n_elements(hole)-1 do begin 
;;;loop over each element of the connected hole                 
                 if hole[k] eq conindex[i,j] then begin ;;;if get a connected match  
;;;want to add all of conindex[i,*] to hole that hasn't already been added
                    for l=0, total(conindex[i,*] gt 0)-1 do begin 
;;;loop over all that the match is connected to                       
                       if total(hole eq conindex[i,l]) eq 0 then begin 
;;;if hasn't already been added to hole                          
                          hole=([hole, conindex[i,l]]) ;;;add to hole                          
                       endif                       
                    endfor                    
                 endif                 
              endfor              
           endfor           
        endif        
     endfor
     
     if n_elements(hole) eq n_elements(zeroindex) then return, 0
  endfor  
  
  ;;;if we get to here then there is a hole as there are unaccounted zeros
  return, 1    
END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; FUNCTION APPLY_SLAVING
;
; Take the reconstructor, remove the rows corresponding to inactive actuators
; and replace them with slaved rows
;
; Inputs: 
;   Hinv: Reconstructor matrix without slaving
;   activeSub: active subaperture map
;   subActuatorMap: maps the actuators adjacent to subapertures
;   actuatorActuatorMap: maps actuators adjacent to actuators
;
FUNCTION APPLY_SLAVING,Hinv,activeSub,subActuatorMap,actuatorActuatorMap  

  whereActiveSub=WHERE(activeSub EQ 1)
  
  activeSubActuatorMap=subActuatorMap[*,whereActiveSub]
  whereNotCorner=WHERE(TOTAL(activeSubActuatorMap,2) EQ 0) ; actuator not touching active subaperture
  
  slaveMatrix=IDENTITY(349.)
  
; check to see if any actuators need to be slaved    
  WHILE whereNotCorner[0] NE -1 DO BEGIN
     neighborActuatorMatrix=actuatorActuatorMap
     neighborActuatorMatrix[whereNotCorner,*]=0.  
     
     FOR c1=0,N_ELEMENTS(whereNotCorner)-1 DO slaveMatrix[*,whereNotCorner[c1]]=neighborActuatorMatrix[*,whereNotCorner[c1]]/(TOTAL(neighborActuatorMatrix[*,whereNotCorner[c1]])+1e-10) ; 1e-10 factor to avoid division by zero
     
     whereNotCorner=WHERE(TOTAL(slaveMatrix,1) EQ 0)   
  ENDWHILE

  Hinv=slaveMatrix##TEMPORARY(Hinv)
  
  RETURN, Hinv
END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
FUNCTION RECONSTRUCT_BAYESIAN,H,weightVector,invcov
  
  weight_matrix=weightVector#REPLICATE(1.,349,1)
  
  alpha=SHOW('ao.RECALPBAY',error=error,status=status,/notrace)
  IF status LT 0 THEN BEGIN
     MESSAGE,/INFO,'PROBLEM READING ao.RECALPBAY'
     MESSAGE,/INFO,error
     alpha=10.
  ENDIF
  
  MESSAGE, /info, 'Alpha coefficient: '+STRING(alpha,format='(f7.3)')
  
  Hinv=invert(TRANSPOSE(H)*weight_matrix##H+alpha*1e-3*invcov+1.)##TRANSPOSE(H)*weight_matrix ; the one is there to penalize piston
  
  RETURN, Hinv
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
FUNCTION RECONSTRUCT_KL,H,weightVector,KL_matrix
  
  weight_matrix=weightVector#REPLICATE(1.,349,1)
  
  KL_modes=SHOW('ao.RECMODKL',error=error,status=status,/notrace)
   IF status LT 0 THEN BEGIN
     MESSAGE,/INFO,'PROBLEM READING ao.RECMODKL'
     MESSAGE,/INFO,error
     KL_modes=250
  ENDIF
  
  Z=KL_matrix[2:KL_modes+1,*] ; modes we want to keep (not tip or tilt)

  AZ=H##Z
  Hinv=Z##invert(transpose(AZ)*weight_matrix##AZ)##transpose(AZ)*weight_matrix

  RETURN, Hinv
END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;         FUNCTION TRIGGER_RECON
;
; This FUNCTION checks whether a new reconstructor needs to be triggered
;
;  
FUNCTION TRIGGER_RECON
  go_recon=0
  
  active_sub_method=SHOW('ao.RECAPSMT',error=error,status=status,/notrace)
  IF STATUS LT 0 THEN BEGIN
     MESSAGE,/INFO,'ERROR READING ao.RECAPSMT'
     MESSAGE,/INFO,error
     active_sub_method=''
  ENDIF  
  
  medianIntensity=SHOW('ao.TLSBMDINT',status=status1, error=err1, /notrace)
  
; if in LGS AO mode, check UTT loop else check DTT loop
  aoopsmode=SHOW('ao.aoopmode',error=error,status=status,/notrace)
  IF aoopsmode EQ 2 THEN BEGIN
     ttlpkwd='ao.utlp' 
     minAllowIntensity=100.     ; minimum allowed intensity to trigger new rec
  ENDIF ELSE BEGIN
     ttlpkwd='ao.dtlp'  
     minAllowIntensity=20.      ; minimum allowed intensity to trigger new rec
  ENDELSE
  
  tiptilt_loop=SHOW(ttlpkwd,status=status2, error=err2, /notrace)
  
  calibrec=(STRUPCASE(active_sub_method) EQ 'ZERO_ANGLE') OR (STRUPCASE(active_sub_method) EQ '304')
  
  IF ((status2 GE 0) AND (status1 GE 0) AND (medianIntensity GE minAllowIntensity) AND tiptilt_loop) OR calibrec THEN BEGIN

;; check if recgo kwd has been set: if it has been set, set go_recon to 1 and clear recgo
     go_recon=SHOW('ao.RECGO',status=status, error=err, /notrace) 
     IF status GE 0 THEN tmp=MODIFY('ao.RECGO',0, status=status, error=err, /notrace)
     ;; check rotator      
     rotmode=SHOW('dcs.rotmode',/notrace)
     CASE rotmode OF
        4: BEGIN ; stationary mode, look at aoprvert
           rotPosActual=SHOW('ao.AOPRVERT',status=statrot, error=errrot, /notrace); in degrees
           rotPosInUse=SHOW('ao.RECANGUSE',status=statrotuse,error=error, /notrace); in degrees   
        END
        ELSE: BEGIN              
           rotPosActual=SHOW('ao.OBRT',status=statrot, error=errrot, /notrace)*180./!pi
           rotPosInUse=SHOW('ao.RECANGUSE',status=statrotuse,error=error, /notrace) ; already in degrees         
        END
     ENDCASE
     diffPos=abs(rotPosActual - rotPosInUse)
     
     IF (diffPos GE 1.) AND (statrot GE 0) AND (statrotuse GE 0) THEN go_recon=1
  ENDIF
  RETURN, go_recon
END

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;   MAIN PROGRAM
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
PRO RECONSTRUCT, data
  DO_NOT_LOAD=0
  
  path2imx='/kroot/rel/ao/qfix/data/ControlParms/SystemMatrix/'  
  path2idlparms='/kroot/rel/ao/qfix/data/IdlParms/'
  
; read in the keywords
  aoopsmode=SHOW('ao.aoopsmode',error=error,status=status,/notrace,/nowait)
  IF aoopsmode EQ 2 THEN LGS_CENTROID_GAIN ; calculate the centroid gain for this pupil angle
  
  active_sub_method=SHOW('ao.RECAPSMT',error=error,status=status,/notrace)
  IF STATUS LT 0 THEN BEGIN
     MESSAGE,/INFO,'ERROR READING ao.RECAPSMT'
     MESSAGE,/INFO,error
     active_sub_method='240'
  ENDIF  
  
  opsmode=SHOW('ao.RECOPSMOD',status=status,error=error,/notrace) 
  IF STATUS LT 0 THEN BEGIN
     MESSAGE,/INFO,'ERROR READING ao.RECOPSMOD'
     MESSAGE,/INFO,error
     opsmode='ON-LINE'
  ENDIF  
  
  IF STRUPCASE(opsmode) EQ 'OFF-LINE' THEN BEGIN
     DO_NOT_LOAD=1
     MESSAGE,/INFO,'Will not load because ao.recopsmod=OFF-LINE'
  ENDIF  

  reconstructor=SHOW('ao.RECMETH',err=err,status=status,/notrace)
  IF STATUS LT 0 THEN BEGIN
     MESSAGE,/INFO,'ERROR READING ao.RECMETH'
     MESSAGE,/INFO,error
     reconstructor='Bayesian'
  ENDIF  
  
; read the plate scale  
  obpsxfs=SHOW('ao.obpsxfs',error=error,status=status,/notrace)
  IF (STATUS LT 0) THEN BEGIN   
     MESSAGE,/INFO,'ERROR READING ao.obpsxfs'
     MESSAGE,/INFO,error
     plat='2.4'
  ENDIF  
  
  CASE obpsxfs OF
     0: plat='2.4'
     1: plat='2.4'
     2: plat='1.0'
     3: plat='0.6'
     ELSE: plat='2.4'
  ENDCASE
  
; load the influence matrix  
  CASE plat OF               
     '1.0': influence_matrix='10.imx'
     '0.6': influence_matrix='06.imx'       
     ELSE: influence_matrix='24.imx'     
  ENDCASE
  
  H=READIMX(path2imx+influence_matrix)
  
; active subaperture methods: 
;   304: all subapertures, even weight
;   1<N<304: N number of subapertures
;   0<N<1.: subapertures with illumination greater than N only
;   Zero_angle: 240 subapertures, even illumination, zero angle  
  
; Measure the subaperture intensity  
  CASE STRUPCASE(active_sub_method) OF
     'ZERO_ANGLE': BEGIN      
        weight=data.weight0
        weightVector=FLTARR(608)
        weightVector[INDGEN(304)*2]=weight[INDGEN(304)]
        weightVector[INDGEN(304)*2+1]=weight[INDGEN(304)]
        activesub=weight
        medianIntensity=1000.
     END
     
     '304': BEGIN
        weight=BYTARR(304)+1    
        weightVector=FLTARR(608)+1.
        activesub=weight
        medianIntensity=1000.
     END
     
     ELSE: BEGIN
        active_sub_method=FLOAT(active_sub_method)
        IF (active_sub_method LE 0) OR (active_sub_method GT 304) THEN BEGIN
           MESSAGE,/INFO,'Active subaperture method '+STRCOMPRESS(active_sub_method,/remove_all)+' is not valid'
           MESSAGE,/INFO,'Using 240 brightest subapertures'
           active_sub_method='240'
        ENDIF
        
        MESSAGE,/info, 'Checking WFS intensity..'
        
        illumination=LONARR(304)
        tiptilt_loop=SHOW('ao.DTLP',status=status2, error=err2, /notrace)
        IF tiptilt_loop EQ 1 THEN illumination=SHOW('ao.tlsbint',error=error,status=status,/notrace)
        medianIntensity=MEDIAN(illumination)
        
        MESSAGE,/info, 'Median intensity ='+STRING(medianIntensity, format='(i5)')   
        
        IF medianIntensity LT 10 THEN BEGIN
           DO_NOT_LOAD=1 
           MESSAGE,/info,'Will not load reconstructor because the median intensity is too low'
        ENDIF 
        
        IF active_sub_method LE 1 THEN BEGIN
           fractional_illum=active_sub_method 
           relative_illum=illumination/medianIntensity
           nActiveSub=TOTAL(relative_illum GT fractional_illum)
           
           IF nActiveSub LT 200 THEN BEGIN 
              nActiveSub=200
              MESSAGE,/INFO,'Fewer than 200 subapertures are active'
              MESSAGE,/INFO,'The number of subapertures is set to 200'
           ENDIF
           
        ENDIF ELSE nActiveSub=active_sub_method        
        
        sortedco=SORT(illumination)           
        activeSub=BYTARR(304)           
        activeSub[sortedco[(304-nActiveSub):303]]=1 ; brightest nActiveSub
        weight=FLTARR(304)
        weight[sortedco[(304-nActiveSub):303]]=illumination[sortedco[(304-nActiveSub):303]]/medianIntensity ; normalize by the median intensity
        weight=weight < 1.
        IF MIN(weight) LT 0 THEN BEGIN
           weight=1.e-3 > weight ; set a floor of 1e-3 (to avoid problems in the LS inversion) and ceiling of 1
           DO_NOT_LOAD=1.
           MESSAGE,/info,'Will not load the reconstrutor because I cannot find '+STRCOMPRESS(nActiveSub,/remove_all)+' subapertures with light' 
        ENDIF
        
        ;;;test to see whether there is a hole in the reconstructor
        IF FINDHOLES(activeSub) EQ 1 THEN BEGIN
           DO_NOT_LOAD=1.
           MESSAGE,/info,'Will not load the reconstrutor because there is a hole in the reconstructor.' 
        ENDIF
        
        weightVector=FLTARR(608)
        
        FOR t1=0,303 DO BEGIN
           weightVector(t1*2.)=(weight(t1)>0.)^1.5
           weightVector(t1*2+1.)=(weight(t1)>0.)^1.5
        ENDFOR
        
        centroid_gain=SHOW('ao.WSCNGN',error=error,status=status,/notrace)
        IF status LT 0 THEN BEGIN
           MESSAGE,/INFO,'Will NOT load the reconstructor because there is an error reading ao.wscngn'
           DO_NOT_LOAD=1
           centroid_gain=FLTARR(608)+1.
        ENDIF
        weightVector=weightVector/centroid_gain^2 ; small cent2arcsec -> small spot -> low noise        
     ENDELSE
  ENDCASE
  
  CASE STRUPCASE(reconstructor) OF
     'BAYESIAN'       : BEGIN
        MESSAGE,/info, 'Calculating Bayesian reconstructor'
        alpha=ADJUST_PARM_B(medianIntensity)
        tmp=MODIFY('ao.RECALPBAY',alpha,error=error,status=status,/notrace)
        Hinv=RECONSTRUCT_BAYESIAN(H,weightVector,data.invcov)  
     END

     'KL'              : BEGIN  ; this is NOT implemented
        MESSAGE,/info, 'Calculating KL reconstructor'
        Hinv=RECONSTRUCT_KL(H,weightVector,data.KL_matrix)
     END

     'SVD'             : BEGIN  ; this is NOT implemented
        MESSAGE,/info, 'Calculating LGS reconstructor'
        Hinv=RECONSTRUCT_SVD(H,activeSub,weightVector)
     END
  ENDCASE
  
; slave 
  slave=SHOW('ao.recslav',status=status,error=error,/notrace)
  CASE slave OF
     -1: BEGIN
        MESSAGE,/INFO,'Could not read ao.recslav keyword'
        MESSAGE,/INFO,'Applying actuator slaving'
        Hinv=APPLY_SLAVING(TEMPORARY(Hinv),activeSub,data.subActuatorMap,data.actuatorActuatorMap) 
        slaving='ON'
     END
     0: BEGIN
        MESSAGE,/INFO,'NOT applying actuator slaving'
        slaving='OFF'
     END
     1: BEGIN
        MESSAGE,/INFO,'Applying actuator slaving'
        Hinv=APPLY_SLAVING(TEMPORARY(Hinv),activeSub,data.subActuatorMap,data.actuatorActuatorMap) 
        slaving='ON'
     END
  ENDCASE
  
; remove piston and tip/tilt
  Hinv=data.ttp##TEMPORARY(Hinv) 
  
; add rows for tip/tilt 
  weight_matrix=weightVector#REPLICATE(1.,2.,1)
  Htilt=FLTARR(2,608)
  Htilt[0,INDGEN(304)*2]=-1.
  Htilt[1,INDGEN(304)*2+1]=-1.
  ttrows=INVERT(TRANSPOSE(Htilt)*weight_matrix##Htilt)##TRANSPOSE(Htilt)*weight_matrix 
  
; remove tip-tilt in centroid space from the DM part of the reconstructor
  msl=IDENTITY(608)-Htilt##ttrows
  Hinv=TEMPORARY(Hinv)##msl
  
; add a row for open loop focus offloading
  focusrow=invert(TRANSPOSE(data.foccents)*weightVector##data.foccents)##TRANSPOSE(data.foccents)*weightVector
  
  Hinvttf=[[-Hinv],[-ttrows],[focusrow]] ; final reconstructor matrix
  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;   Save the new calculated matrix to a file ;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  SPAWN, "date -u | awk '{ print $3$2 }'", root_fn 
  counter=SHOW('ao.RECCYCL',error=error,status=status,/notrace)
  IF STATUS LT 0 THEN counter=0  
  
;; check that there is no NaN in the matrix (could add other quality checks)  
  IF FINITE(MAX(Hinvttf)) EQ 0 THEN BEGIN 
     MESSAGE,/INFO,'Will not load reconstructor because there are NaN values'
     DO_NOT_LOAD=1
  ENDIF 
  
  mr_filename=root_fn+STRING(counter,format="(I4.4)")+'.mr'
  dat_filename=root_fn+STRING(counter,format="(I4.4)")+'.dat'
  
  MESSAGE,/info,'Saving and loading reconstructor!'         
  MESSAGE,/info,'Filename '+mr_filename         
  
  WRITEMR,mr_filename,Hinvttf
  
  prvert=SHOW('ao.aoprvert',error=error,status=status,/notrace)
  alpha=SHOW('ao.recalpbay',error=error,status=status,/notrace)
  
  IF DO_NOT_LOAD EQ 0 THEN BEGIN
     LOADMR,mr_filename,errormsg=errormsg
     IF errormsg NE 'Load Done' THEN BEGIN
        MESSAGE,/INFO,'Failed to load reconstructor'
        MESSAGE,/INFO,errormsg
        SOUNDMESSAGE,'reconFileLoad'
        status=MODIFY('ao.recgo',1,error=error,/notrace)
        RETURN
     ENDIF ELSE BEGIN
; write active subapertures to a keyword     
        tmp=MODIFY('ao.actsubmap',DOUBLE(activesub),error=error,status=status,/notrace)        
        
        rotmode=SHOW('dcs.rotmode',/notrace)
        CASE rotmode OF
           4: rotPosActual=SHOW('ao.AOPRVERT',status=statrot, error=errrot, /notrace)
           ELSE:rotPosActual=SHOW('ao.OBRT',status=statrot, error=errrot, /notrace)*180./!pi
        ENDCASE
        
        IF statrot GE 0 THEN tmp=MODIFY('ao.RECANGUSE',rotPosActual,status=statrotuse,error=error, /notrace)    
        
                                ; need to write a data file with the following info: plate scale, time, active subapertures, regularization parameter or number of modes, active subaperture method, bayesian or modal
        
        OPENW,lun,'/kroot/rel/ao/qfix/data/ControlParms/Recon/'+dat_filename,/get_lun
        PRINTF,lun,STRING('Reconstructor method',format='(A20,''  '')')+STRING(reconstructor,format='(A20)')
        PRINTF,lun,STRING('Alpha',format='(A20,''  '')')+STRING(alpha,format='(A20)')
        PRINTF,lun,STRING('Active subap method',format='(A20,''  '')')+STRING(active_sub_method,format='(A20)')
        PRINTF,lun,STRING('Slaving',format='(A20,''  '')')+STRING(slaving,format='(A20)')
        PRINTF,lun,STRING('Operating mode',format='(A20,''  '')')+STRING(opsmode,format='(A20)')
        PRINTF,lun,STRING('Plate scale',format='(A20,''  '')')+STRING(plat,format='(A20)')
        PRINTF,lun,STRING('Rotator position',format='(A20,''  '')')+STRING(rotPosActual,format='(A20)')
        PRINTF,lun,STRING('Pupil angle',format='(A20,''  '')')+STRING(prvert,format='(A20)')
        FREE_LUN,lun

     ENDELSE
  ENDIF   

END  
  
PRO RECONPROC_EVENT,event

  WIDGET_CONTROL,event.id,get_uvalue=ev
  WIDGET_CONTROL,event.top,get_uvalue=data
    
  IF event.id EQ data.timer_param THEN BEGIN
     go_recon=TRIGGER_RECON()   ; check FOR new recon triggers 
     IF go_recon THEN BEGIN
        counter=SHOW('ao.RECCYCL',status=status,error=error,/notrace)
        IF (STATUS LT 0) OR (counter GT 997) THEN counter=0  
        tmp=MODIFY('ao.RECCYCL',counter+1,error=error,status=status,/notrace) 
        RECONSTRUCT, data
     ENDIF
  ENDIF
  
  WIDGET_CONTROL,event.top,set_uvalue=data     
  WIDGET_CONTROL,event.id,timer=1. ; check and update the screen every second  
  
  path2idlparms='/kroot/rel/ao/qfix/data/IdlParms/'
  SPAWN, 'touch '+path2idlparms+'reconproc.txt'
  
  RETURN
END

PRO RECONPROC  
  path2idlparms='/kroot/rel/ao/qfix/data/IdlParms/'
  
;; find the process ID for the current IDL task  
  currentpid=0L
  SPAWN,'ps',allprocs
  allprocs=allprocs+' '
  idlpidline=WHERE(STRPOS(allprocs,'idl ') NE -1)
  READS,allprocs[idlpidline],currentpid
  OPENW,lun,path2idlparms+'reconproc.txt',/get_lun  
  PRINTF,lun,STRCOMPRESS(currentpid,/remove_all)
  FREE_LUN,lun
  
  zmat=FLTARR(10,608)
  OPENR,lun,'/kroot/rel/ao/qfix/data/ControlParms/zernToCent.dat',get_lun
  READU,lun,zmat
  FREE_LUN,lun
  
  foccents=zmat[3,*]            ;
  
  data={  	$
       mainwid: 1l,   $         ; id of main widget
       timer_param:   0l,   $   ; timer
       foccents:foccents,        $ ; vector of centroids resulting from focus
       invcov: FLTARR(349,349),  $ ; inverse covariance matrix 
       KL_matrix:FLTARR(349,349),$ ; Karhunen-Loeve modes
       ttp: FLTARR(349,349),     $ ; tip-tilt removing matrix
       subActuatorMap: BYTARR(349,304),  $ ; neighboring suabperture matrix 
       actuatorActuatorMap: BYTARR(349,349), $ ; neighboring actuator matrix
       weight0: BYTARR(304)      $ ; zero angle weight matrix
       } 
  
; load the matrix of Karhunen-Loeve modes  
  KL_matrix=FLTARR(349,349)
  OPENR,lun,'/kroot/rel/ao/qfix/data/IdlParms/KL_matrix.mat',/get_lun
  READU, lun, KL_matrix
  FREE_LUN, lun
  data.KL_matrix=KL_matrix
  
; load the inverse covariance matrix
  invcov=FLTARR(349,349)
  OPENR, lun,path2idlparms+'invcov.lsp',/get_lun
  READU,lun,invcov
  FREE_LUN,lun
  data.invcov=invcov
  
; load the tip/tilt removing matrix
  ttp=FLTARR(349,349)
  OPENR,lun,path2idlparms+'ttp.dat',/get_lun
  READU,lun,ttp
  FREE_LUN,lun
  data.ttp=ttp

; load the matrix of whether actuators are touching subapertures 
  subActuatorMap=BYTARR(349,304)
  OPENR, lun,path2idlparms+'actuator_subap.mat',/get_lun
  READU, lun,subActuatorMap
  FREE_LUN,lun
  data.subActuatorMap=subActuatorMap
  
; load the matrix of neighboring actuators  
  actuatorActuatorMap=BYTARR(349,349)
  OPENR, lun,path2idlparms+'neighbour_matrix_actuator.mat',/get_lun
  READU, lun, actuatorActuatorMap
  FREE_LUN,lun
  data.actuatorActuatorMap=actuatorActuatorMap
  
; load the subaperture weights for the zero angle illumination  
  weight0=BYTARR(304)
  OPENR, lun, path2idlparms+'ping0.map', /get_lun ; zero angle map
  READU, lun, weight0
  FREE_LUN, lun     
  data.weight0=weight0
  
; initialize all the keywords    
  tmp=MODIFY('ao.RECOPSMOD', 'MANUAL',err=err,status=status,/notrace)
  tmp=MODIFY('ao.RECMETH', 'BAYESIAN',err=err,status=status,/notrace)
  tmp=MODIFY('ao.RECAPSMT','240',err=err,status=status,/notrace) 
  tmp=MODIFY('ao.RECALPBAY',10,error=error,status=status,/notrace) 
  tmp=MODIFY('ao.recgo',0,error=error,status=status,/notrace)
  tmp=MODIFY('ao.RECSLAV',1,err=err,status=status,/notrace) 
  
  telescope=WHICH_TEL()
  CASE telescope OF
     'Keck I': BEGIN
        tmp=MODIFY('ao.RECMNLGS',1.5,err=err,status=status,/notrace) 
        tmp=MODIFY('ao.RECMXLGS',2.2,err=err,status=status,/notrace) 
     END
     ELSE: BEGIN
        tmp=MODIFY('ao.RECMNLGS',1.5,err=err,status=status,/notrace) 
        tmp=MODIFY('ao.RECMXLGS',3.0,err=err,status=status,/notrace) 
     END
  ENDCASE
  
  data.mainwid=WIDGET_BASE(group_leader=group,column=1,map=0) ; map=0 means widget is hidden
  WIDGET_CONTROL,data.mainwid,/realize
  
  data.timer_param = data.mainwid 
  WIDGET_CONTROL,data.timer_param,timer=1.0
  
; Start the x manager.
  WIDGET_CONTROL,data.mainwid,set_uvalue=data
  XMANAGER,'reconproc',data.mainwid,/no_block
END