restore beam speedup code

fhd_core/beam_modeling/fhd_struct_init_antenna.pro

@@ -34,8 +34,6 @@ zenra=obs.zenra
 zendec=obs.zendec
 obsx=obs.obsx
 obsy=obs.obsy
-;phasera=obs.phasera
-;phasedec=obs.phasedec
 Jdate=obs.Jd0
 IF Keyword_Set(beam_offset_time) THEN Jdate_use=Jdate+beam_offset_time/24./3600. ELSE Jdate_use=Jdate
 frequency_array=(*obs.baseline_info).freq
@@ -123,6 +121,47 @@ valid_i=where(Finite(ra_arr),n_valid)
 ra_use=ra_arr[valid_i]
 dec_use=dec_arr[valid_i]
 
+if ~keyword_set(beam_per_baseline) then undefine, dec_arr, ra_arr
+
+; Split the apply astrometry step into mutiple loops if the image dimension is large to reduce
+; memory footprint
+if keyword_set(conserve_memory) then begin
+  ; calculate bytes required
+  required_bytes = 8.*2*psf_image_dim^2.
+  mem_iter = ceil(required_bytes/mem_thresh)
+endif else mem_iter=1
+if mem_iter GT 1 then begin
+  ; calculate the number of iterations that fit within the array in question
+  mem_iter_array = indgen(mem_iter)
+  mem_iter = max(where((psf_image_dim mod mem_iter_array) EQ 0, n_count))
+  if n_count EQ 0 then mem_iter=1
+  psf_image_dim_use = psf_image_dim / mem_iter
+endif else psf_image_dim_use = psf_image_dim
+
+ra_arr = dblarr(psf_image_dim,psf_image_dim)
+dec_arr = dblarr(psf_image_dim,psf_image_dim)
+
+; Loop through strips in the x dimension to build up the RA and Dec arrays
+for mem_i=0L,mem_iter-1 do begin
+  xvals_celestial=(meshgrid(psf_image_dim_use,psf_image_dim,1)+psf_image_dim_use*mem_i)*psf_scale-psf_image_dim*psf_scale/2.+obsx
+  yvals_celestial=meshgrid(psf_image_dim_use,psf_image_dim,2)*psf_scale-psf_image_dim*psf_scale/2.+obsy
+
+  ;turn off refraction for speed, then make sure it is also turned off in Eq2Hor below
+  apply_astrometry, obs, x_arr=xvals_celestial, y_arr=yvals_celestial, ra_arr=ra_strip, dec_arr=dec_strip, /xy2ad, /ignore_refraction
+  ra_arr[psf_image_dim_use*mem_i:psf_image_dim_use*(mem_i+1)-1,*] = ra_strip
+  dec_arr[psf_image_dim_use*mem_i:psf_image_dim_use*(mem_i+1)-1,*] = dec_strip
+endfor
+
+; Only keep finite pixels (unless full array required later)
+valid_i=where(Finite(ra_arr),n_valid)
+if keyword_set(beam_per_baseline) then begin
+  ra_use=ra_arr[valid_i]
+  dec_use=dec_arr[valid_i]
+endif else begin
+  ra_use=temporary(ra_arr[valid_i])
+  dec_use=temporary(dec_arr[valid_i])
+endelse
+
 ;NOTE: Eq2Hor REQUIRES Jdate_use to have the same number of elements as RA and Dec for precession!!
 ;;NOTE: The NEW Eq2Hor REQUIRES Jdate_use to be a scalar! They created a new bug when they fixed the old one
 Eq2Hor,ra_use,dec_use,Jdate_use,alt_arr1,az_arr1,lat=obs.lat,lon=obs.lon,alt=obs.alt,precess=1,/nutate, refract=0
@@ -136,6 +175,8 @@ antenna.pix_use=ptr_new(pix_use)
 if keyword_set(kernel_window) then begin
   print, 'Applying a modified gridding kernel. Beam is no longer instrumental. Do not use for calibration.'
 
+  ;Get pixels which fall within the horizon
+  horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix)
   xvals_instrument=za_arr*Sin(az_arr*!DtoR)
   yvals_instrument=za_arr*Cos(az_arr*!DtoR)
 
@@ -154,6 +195,9 @@ if keyword_set(kernel_window) then begin
   antenna.pix_window=ptr_new(pix_window)
 endif
 
+if ~keyword_set(pix_use) then horizon_test=where(abs(za_arr) GE 90.,n_horizon_test,complement=pix_use,ncomplement=n_pix)
+antenna.pix_use=ptr_new(pix_use)
+
 ;now, update antenna structure to include gains
 if N_elements(instrument) GT 1 then begin
   for inst_i=0, N_elements(instrument)-1 do begin