; +
; Function to correct MDI magnetogram for cos theta 
;
; P. Gallagher, 1-may-2004
; J McAteer 6-june-2004, made into generic function. 
;   	    	    	 also added in linear / non-linear response.
;J McAteer 20-jan-2005, fixed the cosine correction factor to make it right
; 
; Example: IDL> fits2map,'fd_M_96m_01d.1334.0005.fits', map
;          IDL> cmap = mdi_cos_correct( map )
;
;-

function mdi_cos_correct, map, h_angle=h_angle;, off_centre

    IF KEYWORD_SET(h_angle) THEN h_angle=h_angle ELSE h_angle=60.
    min_correction_factor = cos ((h_angle/90.)*(!pi/2.))
    
    cos_rec   = map.data  ; Corrected array
    ;Remove very large unphysical values, which are limb effects
    ;Also removes large negative value streaks
    IF ((where( cos_rec gt 10000. ))(0) NE -1) THEN cos_rec( where( cos_rec gt 10000. ) ) = 0.
    IF ((where( cos_rec lt -10000.))(0) NE -1) THEN cos_rec( where( cos_rec lt -10000. ) ) = 0.
    sz = size( cos_rec, /dim )
    
    ;Correct for linear and non-liner response of MDI (Berger & Lites 2003)
    flag_nonlin=where(abs(cos_rec) GT 1200.)
    flag_lin=where(abs(cos_rec) LE 1200.)
    
    IF min(flag_nonlin) NE -1 THEN cos_rec(flag_nonlin)=cos_rec(flag_nonlin)*1.9
    IF min(flag_lin) NE -1 THEN cos_rec(flag_lin)=cos_rec(flag_lin)*1.45

    r_sun =  ( pb0r( map.time, /arc, /soho ) )[ 2 ] ; Radius of Sun in arcsec
    
    for i = 0, sz[ 0 ] - 1 do begin 
  
    	for j = 0, sz[ 1 ] - 1 do begin
    
      	iarc = abs(( i - (fix(sz(0)/2.) -1 )) * map.dx + map.xc)  ;calculate x co-ord in arcsec
      	jarc = abs(( j - (fix(sz(1)/2.) -1 )) * map.dy + map.yc)  ;calculate y co-ord in arcsec
      
      	r = sqrt( iarc^2 + jarc^2 )  ;calculate r vector: distance from (0,0)
        r=float(r)
	r_sun=float(r_sun)
	    
      	if ( r / r_sun gt 1. ) then begin    ;account for off-disk pixels
            cos_rec[ i, j ] = 0. 
      	endif else begin
            ;this was peter's old code - incorrect
	    ;cos_theta = cos( [ r / r_sun ] * !pi / 2. )     
            
	    ;this is jma's new code
	    correction_factor = sin (acos(r/r_sun))
	    
	    if correction_factor LT min_correction_factor THEN cos_rec[ i, j ] = randomn(i*j) * 30 ELSE $
	    	cos_rec[ i, j ]= cos_rec[ i, j ] / correction_factor
      	    
	endelse
   
    	endfor
    
    endfor
    
    ;flag if the centre of the region is greater than 60 degrees from disc centre
    ;r = sqrt( map.xc^2 + map.yc^2 )
    ;IF  cos( [ r / r_sun ] * !pi / 2. ) LT 0.5 THEN off_centre= [ r / r_sun ] * 90.
    
    cos_map = map
    add_prop, cos_map, data = cos_rec, /replace
  
    return, cos_map
      
end