apply_astrometry refraction fix (#319)

* Refraction is now a positive, not a double negative

* Changes to `apply_astrometry` refraction for `ad2xy` to those that previously explicitly did not `/ignore_refraction`.
These need to be reviewed

* Changes for `apply_astrometry` `/xy2ad`.

Anything that did not explitcitly set `/ignore_refraction` now has `/refraction` enabled.

Anything that had `/ignore_refraction` now doesn't use the refraction keyword at all.

* Matching what PyFHD does and passes tests

fhd_core/HEALPix/healpix_bin.pro

@@ -57,7 +57,7 @@ IF coord_sys EQ 'galactic' THEN glactc,pix_ra,pix_dec,2000.,pix_ra,pix_dec,2, /d
 IF coord_sys EQ 'equatorial' THEN Hor2Eq,pix_dec,pix_ra,Jdate_use,pix_ra,pix_dec,lat=obs_hpx.lat,lon=obs_hpx.lon,alt=obs_hpx.alt,precess=1,/nutate
 
 ; Transform pixel coordinates to orthoslant coordinates.
-apply_astrometry, obs_hpx, ra_arr=pix_ra, dec_arr=pix_dec, x_arr=xv_hpx, y_arr=yv_hpx, /ad2xy
+apply_astrometry, obs_hpx, ra_arr=pix_ra, dec_arr=pix_dec, x_arr=xv_hpx, y_arr=yv_hpx, /ad2xy, /refraction
 
 ; Select healpix pixels that fall within the field of view.
 hpx_i_use=where((xv_hpx GT 0) AND (xv_hpx LT (dimension_hpx-1)) AND (yv_hpx GT 0) AND (yv_hpx LT (elements_hpx-1)),n_hpx_use) 

fhd_core/HEALPix/healpix_cnv_generate.pro

@@ -23,7 +23,8 @@ IF N_Elements(hpx_radius) EQ 0 THEN BEGIN
         xv_arr=meshgrid(dimension,elements,1)
         yv_arr=meshgrid(dimension,elements,2)
         ;set /ignore_refraction for speed, since we don't need to be exact
-        apply_astrometry, obs, x_arr=xv_arr, y_arr=yv_arr, ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad, /ignore_refraction 
+        ; UPDATE: refraction is ignored by default
+        apply_astrometry, obs, x_arr=xv_arr, y_arr=yv_arr, ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad
         ang_arr=angle_difference(dec_arr,ra_arr,obs.obsdec,obs.obsra,/degree)
         ang_i=where((mask GT 0) AND Finite(ang_arr),n_ang_use)
         IF n_ang_use GT 0 THEN radius=Max(ang_arr[ang_i])

fhd_core/HEALPix/healpix_interpolate.pro

@@ -9,7 +9,7 @@ IF N_Elements(hpx_inds) EQ 0 THEN hpx_inds=Lindgen(n_hpx)
 astr=obs.astr
 dimension=obs.dimension
 elements=obs.elements
-apply_astrometry,obs, x_arr=meshgrid(dimension,elements,1), y_arr=meshgrid(dimension, elements, 2), ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad
+apply_astrometry,obs, x_arr=meshgrid(dimension,elements,1), y_arr=meshgrid(dimension, elements, 2), ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad, /refraction
 radec_i=where(Finite(ra_arr))
 
 IF size(healpix_map,/type) EQ 10 THEN BEGIN ;check if pointer type, and if so allow it to be a pointer array
@@ -33,7 +33,7 @@ vec2ang,pix_coords,pix_dec,pix_ra,/astro
 IF coord_sys EQ 'galactic' THEN glactc,pix_ra,pix_dec,2000.,pix_ra,pix_dec,2, /degree
 IF coord_sys EQ 'equatorial' THEN Hor2Eq,pix_dec,pix_ra,Jdate_use,pix_ra,pix_dec,lat=obs.lat,lon=obs.lon,alt=obs.alt,precess=1,/nutate
 
-apply_astrometry, obs, ra_arr=pix_ra, dec_arr=pix_dec, x_arr=xv_hpx, y_arr=yv_hpx, /ad2xy
+apply_astrometry, obs, ra_arr=pix_ra, dec_arr=pix_dec, x_arr=xv_hpx, y_arr=yv_hpx, /ad2xy, /refraction
 
 hpx_i_use=where((xv_hpx GT 0) AND (xv_hpx LT (dimension-1)) AND (yv_hpx GT 0) AND (yv_hpx LT (elements-1)),n_hpx_use) 
 IF n_hpx_use EQ 0 THEN BEGIN

fhd_core/beam_modeling/fhd_struct_init_antenna.pro

@@ -137,7 +137,8 @@ antenna.psf_scale=psf_scale
 xvals_celestial=meshgrid(psf_image_dim,psf_image_dim,1)*psf_scale-psf_image_dim*psf_scale/2.+obsx
 yvals_celestial=meshgrid(psf_image_dim,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_arr, dec_arr=dec_arr, /xy2ad, /ignore_refraction
+; UPDATE: Refraction is now turned off in apply_astrometry
+apply_astrometry, obs, x_arr=xvals_celestial, y_arr=yvals_celestial, ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad
 undefine, xvals_celestial, yvals_celestial
 valid_i=where(Finite(ra_arr),n_valid)
 ra_use=ra_arr[valid_i]
@@ -168,7 +169,8 @@ for mem_i=0L,mem_iter-1 do begin
   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
+  ; UPDATE: Refraction is now turned off in apply_astrometry
+  apply_astrometry, obs, x_arr=xvals_celestial, y_arr=yvals_celestial, ra_arr=ra_strip, dec_arr=dec_strip, /xy2ad
   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
@@ -200,7 +202,7 @@ if keyword_set(kernel_window) then begin
 
   ;Calculate the phase center in x,y coords
   Eq2Hor,obs.orig_phasera,obs.orig_phasedec,Jdate_use,orig_phasealt,orig_phaseaz,lat=obs.lat,lon=obs.lon,alt=obs.alt
-  apply_astrometry, obs, x_arr=xval_center, y_arr=yval_center, ra_arr=obs.orig_phasera, dec_arr=obs.orig_phasedec, /ad2xy
+  apply_astrometry, obs, x_arr=xval_center, y_arr=yval_center, ra_arr=obs.orig_phasera, dec_arr=obs.orig_phasedec, /ad2xy, /refraction
   xval_center = (90. - orig_phasealt)*sin(orig_phaseaz*!dtor)
   yval_center = (90. - orig_phasealt)*cos(orig_phaseaz*!dtor)
 

fhd_core/calibration/generate_source_cal_list.pro

@@ -134,7 +134,7 @@ IF n_use GT 0 THEN BEGIN
         alpha=spectral_index,extend=catalog.extend)
     endelse
 
-    apply_astrometry, obs, ra_arr=source_list.ra, dec_arr=source_list.dec, x_arr=x_arr, y_arr=y_arr, /ad2xy    
+    apply_astrometry, obs, ra_arr=source_list.ra, dec_arr=source_list.dec, x_arr=x_arr, y_arr=y_arr, /ad2xy, /refraction    
     source_list.x=x_arr
     source_list.y=y_arr
     FOR i=0,7 DO source_list.flux.(i)=catalog.flux.(i)*(freq_use/catalog.freq)^spectral_index
@@ -175,7 +175,7 @@ IF n_use GT 0 THEN BEGIN
         extend_i=where(Ptr_valid(source_list.extend),n_extend)
         FOR ext_i=0L,n_extend-1 DO BEGIN
             extend_list=*source_list[extend_i[ext_i]].extend
-            apply_astrometry, obs, ra_arr=extend_list.ra, dec_arr=extend_list.dec, x_arr=x_arr, y_arr=y_arr, /ad2xy
+            apply_astrometry, obs, ra_arr=extend_list.ra, dec_arr=extend_list.dec, x_arr=x_arr, y_arr=y_arr, /ad2xy, /refraction
             extend_list.x=x_arr
             extend_list.y=y_arr
             *source_list[extend_i[ext_i]].extend=extend_list

fhd_core/deconvolution/fhd_multi_wrap.pro

@@ -51,7 +51,7 @@ FOR j=0L,n_rep-1 DO BEGIN
         normalization=norm_arr2[fi]
         n_pol=obs.n_pol
         astr=obs.astr
-        apply_astrometry, obs, x_arr=meshgrid(obs.dimension,obs.elements,1), y_arr=meshgrid(obs.dimension,obs.elements,2), ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad
+        apply_astrometry, obs, x_arr=meshgrid(obs.dimension,obs.elements,1), y_arr=meshgrid(obs.dimension,obs.elements,2), ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad, /refraction
         dirty_array=Ptrarr(n_pol)
         weights_arr=Reform(weights_arr2[*,fi])
         FOR pol_i=0,n_pol-1 DO BEGIN

fhd_core/deconvolution/source_detection/components2sources.pro

@@ -87,7 +87,7 @@ FOR gi=0L,ng-1 DO BEGIN
     sx=Total(comp_arr[si_g].x*flux_I)/Total(flux_I)
     sy=Total(comp_arr[si_g].y*flux_I)/Total(flux_I)
 
-    apply_astrometry, obs, x_arr=sx, y_arr=sy, ra_arr=sra, dec_arr=sdec, /xy2ad
+    apply_astrometry, obs, x_arr=sx, y_arr=sy, ra_arr=sra, dec_arr=sdec, /xy2ad, /refraction
     source_arr[gi].x=sx
     source_arr[gi].y=sy
     source_arr[gi].ra=sra
@@ -162,7 +162,7 @@ IF Keyword_Set(clean_bias_threshold) THEN BEGIN
             IF n_pix EQ 0 THEN CONTINUE
             sx=Float(pix_i mod dimension)
             sy=Float(Floor(pix_i/dimension))
-            apply_astrometry, obs, x_arr=sx, y_arr=sy, ra_arr=sra, dec_arr=sdec, /xy2ad
+            apply_astrometry, obs, x_arr=sx, y_arr=sy, ra_arr=sra, dec_arr=sdec, /xy2ad, /refraction
             salpha=Replicate(source_arr[si].alpha,n_pix)
             sfreq=Replicate(source_arr[si].freq,n_pix)
             ci=N_Elements(ext_comp)

fhd_core/deconvolution/source_detection/fhd_source_detect.pro

@@ -139,7 +139,7 @@ WHILE n_sources EQ 0 DO BEGIN
         source_box=image_I_flux[sx0-box_radius:sx0+box_radius,sy0-box_radius:sy0+box_radius]
         xcen0=xcen-sx0+box_radius
         ycen0=ycen-sy0+box_radius
-        apply_astrometry, obs, x_arr=xcen, y_arr=ycen, ra_arr=ra, dec_arr=dec, /xy2ad
+        apply_astrometry, obs, x_arr=xcen, y_arr=ycen, ra_arr=ra, dec_arr=dec, /xy2ad, /refraction
         flux_min=Min(source_box[box_radius-local_max_radius:box_radius+local_max_radius,box_radius-local_max_radius:box_radius+local_max_radius])<0
 
         IF flux_interp_flag EQ 0 THEN flux_use=Interpolate(source_box,xcen0,ycen0,cubic=-0.5)>image_I_flux[sx,sy] $

fhd_core/deconvolution/source_detection/fhd_source_detect_healpix.pro

@@ -93,7 +93,7 @@ WHILE n_sources LE 0 DO BEGIN
     si_i_arr=Ptrarr(n_obs)
     FOR obs_i=0L,n_obs-1 DO BEGIN
         si=0L
-        apply_astrometry, obs_arr[obs_i], ra_arr=ra_arr, dec_arr=dec_arr, x_arr=x_arr, y_arr=y_arr, /ad2xy
+        apply_astrometry, obs_arr[obs_i], ra_arr=ra_arr, dec_arr=dec_arr, x_arr=x_arr, y_arr=y_arr, /ad2xy, /refraction
     ;    dist_test=angle_difference(obs_arr[obs_i].obsdec,obs_arr[obs_i].obsra,dec_arr,ra_arr,/degree)
     ;    dist_cut=where(dist_test GT source_alias_radius,n_dist_cut)
 

fhd_core/polarization/fhd_struct_init_jones.pro

@@ -39,16 +39,17 @@ IF N_Elements(mask) EQ 0 THEN mask=Replicate(1.,dimension,elements)
 xvals=meshgrid(dimension,elements,1)
 yvals=meshgrid(dimension,elements,2)
 ;ignore effect of refraction, since we are only determining which pixels to include
-apply_astrometry, obs, ra_arr=ra_arr, dec_arr=dec_arr, x_arr=xvals, y_arr=yvals, /xy2ad, /ignore_refraction
+; UPDATE: refraction is now ignored by default in apply_astrometry
+apply_astrometry, obs, ra_arr=ra_arr, dec_arr=dec_arr, x_arr=xvals, y_arr=yvals, /xy2ad
 inds_use=where(Finite(ra_arr) AND mask,n_pix)
 ra_use=ra_arr[inds_use]
 dec_use=dec_arr[inds_use]
 xv=xvals[inds_use]
 yv=yvals[inds_use]
 
 ;Convert observatory latitude to same coordinate system 
-apply_astrometry, obs, ra_arr=obs.obsra, dec_arr=obs.lat, x_arr=x_t, y_arr=y_t, /ad2xy
-apply_astrometry, obs, dec_arr=lat, x_arr=x_t, y_arr=y_t, /xy2ad, /ignore_refraction
+apply_astrometry, obs, ra_arr=obs.obsra, dec_arr=obs.lat, x_arr=x_t, y_arr=y_t, /ad2xy, /refraction
+apply_astrometry, obs, dec_arr=lat, x_arr=x_t, y_arr=y_t, /xy2ad
 
 obs_temp = fhd_struct_update_obs(obs, beam_nfreq_avg=obs.n_freq) ; Use only one average Jones matrix, not one per frequency
 antenna=fhd_struct_init_antenna(obs_temp,beam_model_version=beam_model_version,psf_resolution=1.,$

fhd_core/polarization/rotate_jones_matrix.pro

@@ -15,7 +15,7 @@ endfor
 ;Generate a grid of RA and Dec coordinates matching the image
 xvals=meshgrid(obs.dimension,obs.elements,1)
 yvals=meshgrid(obs.dimension,obs.elements,2)
-apply_astrometry, obs, ra_arr=ra_arr, dec_arr=dec_arr, x_arr=xvals, y_arr=yvals, /xy2ad
+apply_astrometry, obs, ra_arr=ra_arr, dec_arr=dec_arr, x_arr=xvals, y_arr=yvals, /xy2ad, /refraction
 ; Pixels beyond the horizon will have NAN RA and Dec values. Skip them.
 inds_use=where(Finite(ra_arr),n_pix)
 ra_use=ra_arr[inds_use]

fhd_core/source_modeling/fhd_galaxy_model.pro

@@ -13,7 +13,7 @@ elements=obs.elements
 astr=obs.astr
 degpix=obs.degpix
 n_pol=obs.n_pol
-apply_astrometry, obs, x_arr=meshgrid(dimension,elements,1), y_arr=meshgrid(dimension,elements,2), ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad
+apply_astrometry, obs, x_arr=meshgrid(dimension,elements,1), y_arr=meshgrid(dimension,elements,2), ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad, /refraction
 IF N_Elements(jones) EQ 0 THEN jones=fhd_struct_init_jones(obs,file_path_fhd=file_path_fhd,/restore)
 
 ;i_use=where(Finite(ra_arr))

fhd_core/source_modeling/source_dft_multi.pro

@@ -72,7 +72,7 @@ IF keyword_set(gaussian_source_models) then begin
 
         ;Calculate the x,y pixel coordinates of the gaussian sources
         apply_astrometry, obs, x_arr=gaussian_x_vals, y_arr=gaussian_y_vals, ra_arr=[[obs.obsra-0.5*gaussian_ra], $
-          [obs.obsra+0.5*gaussian_ra]], dec_arr=[[obs.obsdec-0.5*gaussian_dec],[obs.obsdec+0.5*gaussian_dec]], /ad2xy
+          [obs.obsra+0.5*gaussian_ra]], dec_arr=[[obs.obsdec-0.5*gaussian_dec],[obs.obsdec+0.5*gaussian_dec]], /ad2xy, /refraction
         gaussian_x = (gaussian_x_vals[*,1]-gaussian_x_vals[*,0]) / gaussian_ra_corr
         gaussian_y = (gaussian_y_vals[*,1]-gaussian_y_vals[*,0]) / gaussian_dec_corr
 

fhd_core/visibility_manipulation/vis_delay_filter.pro

@@ -40,7 +40,7 @@ pro vis_delay_filter, vis_model_arr,  params, obs
 
   ; Phase to zenith -- easier calculations of the location of the horizon when phased to zenith
   dimension=obs.dimension
-  apply_astrometry,obs,ra_arr=obs.zenra, dec_arr=obs.zendec, x_arr=x_use, y_arr=y_use, /ad2xy
+  apply_astrometry,obs,ra_arr=obs.zenra, dec_arr=obs.zendec, x_arr=x_use, y_arr=y_use, /ad2xy, /refraction
   dx=obs.obsx - obs.zenx
   dy=obs.obsy - obs.zeny
   dx*=(2.*!Pi/dimension)

fhd_output/fhd_quickview.pro

@@ -113,7 +113,8 @@ astr_out=obs_out.astr
 horizon_mask=fltarr(dimension,elements)+1.
 ;IF Keyword_Set(image_mask_horizon) THEN BEGIN
     ;set /ignore_refraction for speed since we're just finding pixels to mask
-    apply_astrometry, obs_out, x_arr=meshgrid(dimension,elements,1), y_arr=meshgrid(dimension,elements,2), ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad, /ignore_refraction
+    ; UPDATE: Refraction is ignored by default in apply astrometry
+    apply_astrometry, obs_out, x_arr=meshgrid(dimension,elements,1), y_arr=meshgrid(dimension,elements,2), ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad
     horizon_test=where(Finite(ra_arr,/nan),n_horizon_mask)
     IF n_horizon_mask GT 0 THEN horizon_mask[horizon_test]=0
 ;ENDIF
@@ -199,14 +200,14 @@ IF model_flag THEN IF skymodel.n_sources GT 0 THEN BEGIN
     source_array=skymodel.source_list
     source_arr_out=source_array
 
-    apply_astrometry, obs_out, ra_arr=source_array.ra, dec_arr=source_array.dec, x_arr=sx, y_arr=sy, /ad2xy
+    apply_astrometry, obs_out, ra_arr=source_array.ra, dec_arr=source_array.dec, x_arr=sx, y_arr=sy, /ad2xy, /refraction
     source_arr_out.x=sx & source_arr_out.y=sy
 
     extend_test=where(Ptr_valid(source_arr_out.extend),n_extend)
     IF n_extend GT 0 THEN BEGIN
         FOR ext_i=0L,n_extend-1 DO BEGIN
             component_array_out=*source_array[extend_test[ext_i]].extend
-            apply_astrometry, obs_out, ra_arr=component_array_out.ra, dec_arr=component_array_out.dec, x_arr=cx, y_arr=cy, /ad2xy
+            apply_astrometry, obs_out, ra_arr=component_array_out.ra, dec_arr=component_array_out.dec, x_arr=cx, y_arr=cy, /ad2xy, /refraction
             component_array_out.x=cx & component_array_out.y=cy
 
             IF Total(component_array_out.flux.(0)) EQ 0 THEN BEGIN
@@ -229,14 +230,14 @@ IF size(cal, /type) EQ 8 THEN IF cal.skymodel.n_sources GT 0 THEN BEGIN
     cal_source_array = cal.skymodel.source_list
     cal_source_arr_out=cal_source_array
 
-    apply_astrometry, obs_out, ra_arr=cal_source_array.ra, dec_arr=cal_source_array.dec, x_arr=sx, y_arr=sy, /ad2xy
+    apply_astrometry, obs_out, ra_arr=cal_source_array.ra, dec_arr=cal_source_array.dec, x_arr=sx, y_arr=sy, /ad2xy, /refraction
     cal_source_arr_out.x=sx & cal_source_arr_out.y=sy
 
     extend_test=where(Ptr_valid(cal_source_arr_out.extend),n_extend)
     IF n_extend GT 0 THEN BEGIN
         FOR ext_i=0L,n_extend-1 DO BEGIN
             component_array_out=*cal_source_array[extend_test[ext_i]].extend
-            apply_astrometry, obs_out, ra_arr=component_array_out.ra, dec_arr=component_array_out.dec, x_arr=cx, y_arr=cy, /ad2xy
+            apply_astrometry, obs_out, ra_arr=component_array_out.ra, dec_arr=component_array_out.dec, x_arr=cx, y_arr=cy, /ad2xy, /refraction
             component_array_out.x=cx & component_array_out.y=cy
 
             cal_source_arr_out[extend_test[ext_i]].extend=Ptr_new(/allocate)

fhd_utils/format_conversion/globalskymodel_read.pro

@@ -10,7 +10,7 @@ elements=obs.elements
 astr=obs.astr
 degpix=obs.degpix
 n_pol=obs.n_pol
-apply_astrometry, obs, x_arr=meshgrid(dimension, elements, 1), y_arr=meshgrid(dimension, elements, 2), ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad
+apply_astrometry, obs, x_arr=meshgrid(dimension, elements, 1), y_arr=meshgrid(dimension, elements, 2), ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad, /refraction
 
 file_path_base=filepath('',root=rootdir('FHD'),sub='catalog_data')
 IF Keyword_Set(haslam_filtered) THEN BEGIN
@@ -45,7 +45,7 @@ IF Keyword_Set(haslam_filtered) THEN BEGIN
     GlactC,pix_ra,pix_dec,2000.,pix_gl,pix_gb,2,/degree
     Eq2hor,pix_ra,pix_dec,obs.JD0,pix_alt,pix_az,lat=obs.lat,lon=obs.lon,alt=obs.alt, /precess, /nutate, refract=0
 
-    apply_astrometry,obs, ra_arr=pix_ra, dec_arr=pix_dec, x_arr=xv_hpx, y_arr=yv_hpx, /ad2xy
+    apply_astrometry,obs, ra_arr=pix_ra, dec_arr=pix_dec, x_arr=xv_hpx, y_arr=yv_hpx, /ad2xy, /refraction
     hpx_i_use=where((xv_hpx GT 0) AND (xv_hpx LT (dimension-1)) AND (yv_hpx GT 0) AND (yv_hpx LT (elements-1)) AND (pix_alt GT 0),n_hpx_use,complement=hpx_i_cut) 
 ;    hpx_i_use=hpx_i_cut
     IF n_hpx_use EQ 0 THEN RETURN,Ptrarr(n_pol)

fhd_utils/format_conversion/mrc_catalog_read.pro

@@ -47,7 +47,7 @@ IF Keyword_Set(obs) THEN BEGIN
     angs=angle_difference(dec0,ra0,dec,ra,/degree)
     i_use=where(Abs(angs) LE 60.,n_use)
     IF n_use GT 0 THEN BEGIN
-        apply_astrometry, obs, ra_arr=ra[i_use], dec_arr=dec[i_use], x_arr=x_arr, y_arr=y_arr, /ad2xy
+        apply_astrometry, obs, ra_arr=ra[i_use], dec_arr=dec[i_use], x_arr=x_arr, y_arr=y_arr, /ad2xy, /refraction
         mrc_cat[i_use].x=x_arr
         mrc_cat[i_use].y=y_arr
     ENDIF

fhd_utils/horizon_mask.pro

@@ -6,8 +6,9 @@ elements = obs.elements
 horizon_mask=intarr(dimension,elements)+1
 
 ;ignore refraction since it's only being used to calculate a pixel mask for the horizon
+; UPDATE: refraction is now ignored by default in apply_astrometry 
 apply_astrometry, obs, x_arr=meshgrid(dimension,elements,1), y_arr=meshgrid(dimension,elements,2), $
-    ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad, /ignore_refraction
+    ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad
 nan_i=where(Finite(ra_arr,/nan),n_nan,complement=horizon_i)
 IF n_nan GT 0 THEN horizon_mask[nan_i]=0
 ra_use=ra_arr[horizon_i]

fhd_utils/modified_astro/apply_astrometry.pro

@@ -21,14 +21,14 @@
 ;    
 ;    ad2xy - set to convert from celestial to pixel coordinates
 ;    
-;    ignore_refraction - set to ignore refraction (for backwards compatibility and testing)
+;    refraction - set to take refraction into account (it's false by default for speed, backwards compatibility and testing)
 ;    
 ;    _Extra - passes any additional keywords (like temperature (Kelvin) or pressure (millibars) to co_refract
 ;
 ; :Author: Ian
 ;-
 PRO apply_astrometry,obs, x_arr=x_arr, y_arr=y_arr, ra_arr=ra_arr, dec_arr=dec_arr, xy2ad=xy2ad, ad2xy=ad2xy, $
-    astr=astr, JDate=JDate, ignore_refraction=ignore_refraction, _Extra=extra
+    astr=astr, JDate=JDate, refraction=refraction, _Extra=extra
 
 IF not (Keyword_Set(xy2ad) OR Keyword_Set(ad2xy)) THEN $
     message, "ERROR! At least one direction keyword (xy2ad or ad2xy) must be set."
@@ -39,13 +39,16 @@ IF (Keyword_Set(xy2ad) AND Keyword_Set(ad2xy)) THEN $
 IF size(obs,/type) EQ 8 THEN BEGIN 
     astr=obs.astr
     JDate=obs.JD0
-ENDIF ELSE ignore_refraction=1
+ENDIF ELSE refraction=0
+
+if Keyword_Set(refraction) THEN refraction=1 ELSE refraction=0
+
 precess_forward=0
 precess_reverse=0
 IF Keyword_Set(xy2ad) THEN BEGIN
     xy2ad, x_arr, y_arr, astr, ra_arr, dec_arr
 
-    IF ~Keyword_Set(ignore_refraction) THEN BEGIN   
+    IF Keyword_Set(refraction) THEN BEGIN   
         i_nan=where(Finite(ra_arr,/nan),n_nan,complement=i_use)
         IF n_nan EQ 0 THEN BEGIN
             Eq2Hor,ra_arr, dec_arr, JDate, alt_arr, az_arr, nutate=precess_reverse,precess=precess_reverse,aberration=precess_reverse, refract=0, lon=obs.lon, alt=obs.alt, lat=obs.lat
@@ -78,7 +81,7 @@ IF Keyword_Set(ad2xy) THEN BEGIN
         dec_arr_new[horizon_i] = !Values.F_NAN
         i_nan=where(Finite(ra_arr_new,/nan),n_nan,complement=i_use)
     ENDIF
-    IF ~Keyword_Set(ignore_refraction) THEN BEGIN    
+    IF Keyword_Set(refraction) THEN BEGIN    
         ra_arr_new = ra_arr
         dec_arr_new = dec_arr
         IF n_nan EQ 0 THEN BEGIN

fhd_utils/modified_astro/pixel_area.pro

@@ -5,7 +5,8 @@ elements=obs.elements
 xvals=meshgrid(dimension,elements,1)
 yvals=meshgrid(dimension,elements,2)
 ;refraction has negligible effect on the size of pixels, so ignore to save time
-apply_astrometry, obs, x_arr=xvals, y_arr=yvals, ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad, /ignore_refraction
+; UPDATE: refraction is now ignored by default in apply_astrometry
+apply_astrometry, obs, x_arr=xvals, y_arr=yvals, ra_arr=ra_arr, dec_arr=dec_arr, /xy2ad
 ra_i_nan=where(Finite(ra_arr,/nan),n_nan,complement=ra_i_use)
 ra_vals=ra_arr[ra_i_use]
 dec_vals=dec_arr[ra_i_use]