;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;
function vignettage,xc,yc
;;  created by M.B on 02/02/94 at LAS
;;  modified by M.B on 04/03/94 at LAS
;;  construction of a normalized model of vignetting function for LASCO
;;  This model can be used for vignetting of solar corona frames
;;  inputs :
;;
;;          xc,yc   center coordinates of internal occultor
;;
;;  outputs :
;;          model         name of the resultant image
;;
;;  restrictions :
;;          the generated model is 512*512 pixels
;
;-------------------------------------------------------------------------------
;
scale1 = 2.*0.021*864./364.
scale2 = 2.*0.021*820./364.
RD2 = (31.94*1.08/2.)/scale2
RA3i = (21.*0.96/2.)/scale2
RA3e = (21.*0.96/2.)/scale1
RA2 = 7.19*(32.4/(7.52/1.08))/scale1
print,'>>> begin construction of the equivalent occulting disk D2 <<<'
GENER_MASK,D2MASK,512,xc,yc,RD2,1.,0.
print,'>>> begin construction of the equivalent pylon of D2 <<<'
PYLON = BAND_MASK(512,79)
PYLON(214:297,256:511) = 0.
PYLON = rot_int(PYLON,47.,1,255.5,255.5)
D2MASK = D2MASK + PYLON
D2MASK = D2MASK < 1.000001
D2MASK = 1.-D2MASK
print,'>>> begin construction of the equivalent entrance pupil A3 <<<'
GENER_MASK,A3MASKI,512,xc,yc,RA3i,1.,0.
GENER_MASK,A3MASKE,512,xc,yc,RA3e,1.,0.
print,'>>> begin construction of the equivalent internal diaphragm A2 <<<'
GENER_MASK,A2MASK,700,349.5,349.5,RA2,1.,0.
A2MASK = EXTRAC(A2MASK,94,94,512,512)
dx = 255.5 - xc
dy = 255.5 - yc
x1 = 255.5 + dx
y1 = 255.5 + dy
A2MASK = ROT_INT(A2MASK,0.,1.,x1,y1)
;
print,'>>> begin construction of internal diaphragm <<<'
print,'>>> begin convolution between entrance pupil & occulting system <<<'
diaph_int = n_elements(A3MASKI)*float(dfti(dft(D2MASK)*dft(A3MASKI)))
diaph_int = diaph_int/max(diaph_int)
;
print,'>>> begin construction of external diaphragm <<<'
print,'>>> begin convolution between entrance pupil & internal diaphragm <<<'
diaph_ext = n_elements(A3MASKE)*float(dfti(dft(A2MASK)*dft(A3MASKE)))
diaph_ext = diaph_ext/max(diaph_ext)
;
model = diaph_int*diaph_ext
model = model > 0.
;
print,'>>> statistics & visualization of the resulting image <<<'
;
stat_ima,model,0,0,511,511
wset,0
loadct,3
tvscl,bytscl(model,min=0.,max=max(model))
wset,1
plot,model(*,yc),yrange=[0.,max(model)],linestyle=0
oplot,model(xc,*),linestyle=1
;
print,'end of job'
;
return,model
end
;
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; NAME:                        crevig.pro
; PURPOSE:                     construction of the vignetting function for
;                              LASCO/C2 coronagraph using the bidimensionnal
;                              convolution filtering method
; CATEGORY:                    Processing high level
; CALLING SEQUENCE:            crevig,ima_name,off_x,off_y
; INPUTS:                      ima_name      name of the simulation result
;                              off_x,off_y   off-centering of the occultor
;
; OPTIONAL INPUT PARAMETERS:   None
; KEYWORD PARAMETERS:          None
; OUTPUTS:                     image of given name
;
; OPTIONAL OUTPUT PARAMETERS:  None
; COMMON BLOCKS:               None
; SIDE EFFECTS:                None
; RESTRICTIONS:                This procedure is adapted to 512*512 frames
; PROCEDURE:
; MODIFICATION HISTORY:   defined by M.B 02/07/94
;
; SCCS variables for IDL use
;
; @(#)crevig.pro  v 1.0 02/07/94 :LAS
;
;-------------------------------------------------------------------------------
;
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
function lunule,R1,R2,d
;-------------------------------------------------------------------------------
;
if (R1 gt R2) then begin $
  A = R1
  B = R2
endif else begin $
  A = R2
  B = R1
endelse
;
r = abs(d)
A2 = A*A
B2 = B*B
SEP = 1.
;
if (r gt A-B) then begin $
  if (r ge A+B) then begin $
    SEP = 1.
  endif else begin $
    eta = (A2-B2)/(2.*r)
    zeta = sqrt(.5*(A2+B2)-(eta^2 + (.5*r)^2))
    tetha_A = acos((.5*r+eta)/A)
    seg_A = A2*tetha_A-(eta+.5*r)*zeta
    if (eta gt .5*r) then begin $
      tetha_B = acos((eta-.5*r)/B)
      seg_B = B2*tetha_B-(eta-.5*r)*zeta
    endif else begin $
      tetha_B = 3.14159263-acos((.5*r-eta)/B)
      seg_B = B2*tetha_B-(eta-.5*r)*zeta
    endelse
    area = 3.14159263*B2
    SEP = (seg_B-seg_A)/area
  endelse
endif else begin $
  SEP = 0.
endelse
;
return,SEP
end
;
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;
pro crevig,ima_name,off_x,off_y
;
;-------------------------------------------------------------------------------
;
RA11 = 10.1
RD01 = 17.245*0.992
RA20 = 35.388
diaph_int = fltarr(512,512)
diaph_ext = fltarr(512,512)
ima_name = fltarr(512,512)
;
step_int = 2.*0.021*820./364.
step_ext = 2.*0.021*864./364.
cx_int = off_x*step_int
cy_int = off_y*step_int
start_int = -step_int*255.5
cx_ext = off_x*step_ext
cy_ext = off_y*step_ext
start_ext = -step_ext*255.5
;
for i=1,512 do begin $
  y = (i-1)*step_int + start_int
  for j=1,512 do begin $
  x = (j-1)*step_int + start_int
    D2 = (x-cx_int)^2 + (y-cy_int)^2
    if (D2 gt 0.) then begin $
      D1 = sqrt(D2)
    endif else begin $
      D1 = 0.
    endelse
    diaph_int(i-1,j-1) = lunule(RD01,RA11,D1)
  endfor
endfor
;
for i=1,512 do begin $
  y = (i-1)*step_ext + start_ext
  for j=1,512 do begin $
  x = (j-1)*step_ext + start_ext
    D2 = (x-cx_ext)^2 + (y-cy_ext)^2
    if (D2 gt 0.) then begin $
      D1 = sqrt(D2)
    endif else begin $
      D1 = 0.
    endelse
    diaph_ext(i-1,j-1) = lunule(RA20,RA11,D1)
  endfor
endfor
;
ima_name = diaph_int*(1.-diaph_ext)
ima_name = ima_name > 0.
;
stat_ima,ima_name,0,0,511,511
wset,0
loadct,3
tvscl,bytscl(ima_name,min=0.,max=1.)
wset,1
plot,ima_name(*,255),linestyle=0
oplot,ima_name(255,*),linestyle=1
;
return
end