;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; NAME:                        lux_diff.pro
; PURPOSE:                     make an image of the scattering function
; CATEGORY:                    Processing high level
; CALLING SEQUENCE:            lux_diff,ima_name,Xc,Yc,OPTION,A,B
; INPUTS:                      ima_name                       Name of result
; OPTIONAL INPUT PARAMETERS:   Xc,Yc                          center's coords
;                              OPTION                         function's chars
; KEYWORD PARAMETERS:          None
; OUTPUTS:                     ima_name                       scattering funct
; OPTIONAL OUTPUT PARAMETERS:  None
; COMMON BLOCKS:               None
; SIDE EFFECTS:                None
; RESTRICTIONS:                Applications limited to 512*512 frames
; PROCEDURE:
; MODIFICATION HISTORY:   defined by M.B 10/11/93
;
; SCCS variables for IDL use
;
; @(#)lux_diff.pro  1.0 10/11/93 :LAS
;
;----------------------------------------------------------------------------
;
pro lux_diff,ima_name,Xc,Yc,OPTION,A,B
;
img = MAKE_ARRAY(512,512)
;
for x=0,511 do begin $
  delta_x = ABS(x-255.)
  deltax2 = delta_x * delta_x
for y=0,511 do begin $
  delta_y = ABS(y-255.)
  deltay2 = delta_Y * delta_Y
  xr = SQRT(ABS(deltax2 + deltay2))
  xr = xr/40.73
  IF (xr GT 1.) THEN BEGIN $
  case OPTION of
  1:      img(x,y) = A*(xr)^(-B)
  2:      img(x,y) = 0.53703*(xr)^(-1)
  3:      img(x,y) = 1.58E03*(xr)^(-1.37)   ; DURST's scattering function (1982)
  4:      img(x,y) = 0.53703*(xr)^(-2)
  5:      OPTION = 5
  ENDCASE
  ENDIF
  endfor
endfor
  IF (OPTION EQ 5) THEN BEGIN $; KOUTCHMY-KOUTCHMY 's scattering function (1974)
   for x=0,511 do begin $
       delta_x = ABS(x-255.)
       deltax2 = delta_x * delta_x
     for y=0,511 do begin $
       delta_y = ABS(y-255.)
       deltay2 = delta_Y * delta_Y
       xr = SQRT(ABS(deltax2 + deltay2))
       xr = xr/40.73
       IF (xr LE 2.5 AND xr GT 1.) THEN BEGIN $
       img(x,y) = 10.^(-2.2325-(4.664284/xr)+(1.28571/(xr^2)))
       ENDIF
       IF (xr GT 2.5) THEN BEGIN $
       img(x,y) = ((2.5/xr)^3)*(10.^(-2.2325-(4.664284/2.5)+(1.28571/(2.5^2))))
       ENDIF
      endfor
    endfor
  ENDIF
imab = img
Dx = FIX(Xc-255.)
Dy = FIX(Yc-255.)
imab = SHIFT(imab,Dx,Dy)
ima_name = imab
;
MX = MAX(ima_name)
MN = MEDIAN(ima_name)
MN = MN/10.
mu = make_array(2,512)
mu(0,*) = (indgen(512)-Xc)/40.73
mu(1,*) = ima_name(*,Yc)
WSET,1
ERASE
LOADCT,12
TVSCL,BYTSCL(ima_name,MIN=MN,MAX=MX)
WSET,0
ERASE
PLOT_IO,mu(0,*),mu(1,*),YRANGE=[MN,MX],LINESTYLE=0
STAT_IMA,ima_name,0,0,511,511
RETURN
END