pro cal_slope,image,mn,mx,angle
;
;+
;*****************************************************
;**NAME: CAL_SLOPE
;**CATEGORY: IRAS Image Processing routine: calculates the slope of 
; 	     periodic stripes in IRAS skyflux images.
;**CALLING SEQUENCE:  cal_slope, image, mn, mx, angle
;**PARAMETERS:  IMAGE, MN, MX, ANGLE
;	IMAGE:  (REQ) (I) (RS) (0,2.. either string or 2-D array)
;		A median filtered version of the background 
;		(this is the input)
;	MN, MX: (OPT) (O) (R) (0)
;		if present, will define the minimum and maximum value
;		of the array 'angle' (respectively).  The angles are
;		defined as positive if pointing north/east, as negative
;		if pointing north/west, and as zero, if straight up
;		(vertical).  If both parameters are absent, angles
;		are estimated by the program.
;	ANGLE: (REQ) (O) (R) (1)
;		This is the array which contains the angles that the
;		stripes make with the y-axis.
;**EXAMPLES:
;	is called within IRAS_DESTRIPE
;**INTERACTIVE INPUT:
;	none
;**COMMON BLOCKS:
;	none
;**SUBROUTINES CALLED:
;	none
;**FILES USED:
;	none
;**RESTRICTIONS:
;	Should work for versions 1 and 2 of IDL.
;**SUBROUTINES CALLED: 
;	FRQ
;**PROCEDURE:
;	The procedure determines the slope of the stripes in IRAS 
;	skyflux images, defining positive angles as pointing north-
;	west and negative angles north-east.
;	The stripes are calculated as one sigma away from the median va-
;	lue of the array 'angle'.
;**MODIFICATION HISTORY:
;	July 1990, Written by Vasudevan Subramanian
;	September 1990, Header written by Jose Ruiz
;
;********************************************
;-
;
;STEP 1: initializations and a call to the period estimation program

frq,image,period
vert_space = nint(period/1.5)
shift_lim  = nint(0.4*period) 
                                 ; In the ideal case make this just
                                 ; less than period/2
err   = fltarr(2*shift_lim+1)
errprime = fltarr(2*shift_lim+1)
angle = fltarr(25)
loc_minval  = fltarr(150)
loc_minindex = lonarr(150)
row1  = lonarr(5)
row2  = lonarr(5)
col1  = lonarr(5)
index = 0
noisy = 0                        ; Counter for number of cases
                                 ; where no minimum is yielded
;STEP 2: compute the sums of absolute errors
;        for each sector of the image

for l = 0,4 do begin
 col1(l) = 50*l + 100            ; column locations for error computations
                                
 for m = 0,4 do begin 
  row1(m) = 50*m + 100          ; row locations for error computations
  x_shift = 0                  
                        
   for k =0,1 do begin           ;test for +ve(k=1) as well as 
                                 ;-ve(k=0) slopes

    
    row1(m) = row1(m) + k*vert_space
    row2(m) = row1(m) + (1-(2*k))*vert_space

    for i=k*(shift_lim+1),(k+1)*shift_lim do begin
    
     errprime(i) = 0.0                  ;sum up the errors
     dummy_var   = i - k*shift_lim      ;for various shifts          
   
     for j=col1(l),col1(l)+nint(1.5*period)  do begin
      errprime(i) = errprime(i) + $
                    abs(image(j,row2(m))-image(dummy_var+j,row1(m)))
     endfor                      ; (j)
 
    endfor                       ; (i)

   endfor                        ; (k)
 
;STEP 3: find the minimum value in the error sequence
;        and the corresponding index, simultaneously
;        determining the direction (+ve) or (-ve)

;STEP 3.1: restring the error vector

   for i=0,shift_lim do begin
    err(i+shift_lim)=errprime(i)
   endfor
   for i=0,shift_lim-1 do begin
    err(i) = errprime(2*shift_lim-i)
   endfor

   ;plot,err

;STEP 3.2: locate the local minima and then the global minimum in 
;          the error array. When none are found, a counter named 
;          noisy is incremented.
;          This eliminates the possibility of incorrect estimates,
;          since obtaining a minimum at the end points; err(0),
;          or  err(2*shift_lim) indicates that the true minimum 
;          has not been found.

   g=0
   loc_minval(0)   = 0.0
   loc_minindex(0) = 0
 

   for q = 1,2*shift_lim-1 do begin  ; obtain the local minima
    if ((err(q) le err(q+1)) and (err(q) le err(q-1))) then begin
     loc_minval(g) = err(q)
     loc_minindex(g) = q
     g = g+1
    endif   
   endfor                  ; (q)

   glob_min = loc_minval(0)
   x_shift  = loc_minindex(0)

   if ( g gt 1) then begin
    for c = 1,g-1 do begin
     if (glob_min ge loc_minval(c)) then begin
      glob_min = loc_minval(c)
      x_shift  = loc_minindex(c)
     endif 
    endfor                 ; (c)
   endif
   
   minerr = min(err)
   if ( (g eq 0) or (glob_min gt minerr) ) then begin
    ;print,'No minimum found'
   endif else begin
                           ; compute the slopes with the signs
    x_shift = shift_lim - x_shift 
    ;print,x_shift
    y_shift = abs(row1(m) - row2(m))
    argument = float(x_shift)/float(y_shift)   
    angle(index) = atan(argument)
    angle(index) = angle(index)*!radeg
  
    ;print,argument
    ;print,angle(index)
    index = index + 1
   endelse
   
   
 endfor         ; (m) end of loops for angle computations
endfor          ; (l) at various points

angles=angle(0:index-2) ;"angles" contains all reasonable numbers in "angle"
sortangles=angle(sort(angles)) ;sort angles in ascending order
stand=stdev(angles) ;standard deviation of array "angles"
mx=sortangles((index-1)/2)+stand 
mn=sortangles((index-1)/2)-stand
return
end