pro pol_std_analysis
;
;   analyze polarization for standard (defocussed) star
;
;   performs analysis via three routes:
;   1. by forming the Qi and Ui quantities, for which there are normally 8 each in each observation,
;     then finding the means and mean uncertainties for these quantities for each aperture.
;     Finally, the aperature with the minimum uncertainty in mean Q has this Q value selected,
;     and the same for mean U (though the apertures need not be the same). Polarization, PA, and
;     Ricean corrected P and sPA are computed. From analysis of 22 E22 observations, this approach
;     best matches the variation of P to its sP and variation of PA to its sPA. As of 1/1/2006,
;     this is the source of the output values.
;   2. by folding the Fluxes onto each other for the fluxes taken at the "same" HWP phase (there are
;     four full cycles of this phase per HWP revolution). In this analysis, the mean flux is computed
;     from the photometric error-weighted values and the mean error of the flux is computed from
;     the dispersion of the flux values reduced by sqrt N-1 of the samples. This analysis either
;     returns uncertainties that are too small (if 1/sqrt(sum of weights) is used) or too large
;     (if dispersion/sqrt(N-1) is used)), probably due to improper or excessive accounting of sky
;     drifts and noise.
;   3. by sinusoidal fitting of the folded fluxes from (2.) above. When the uncertainties are very
;     small, all three methods recover identical results. However, as uncertainties grow, the sinusoidal
;     constants are not well constrained and they wander. Also, for low polarizations, the sinusoidal
;     fitting will always find higher polarizations than are real.
;   As of 1/1/2006, I have moved to save only the results from (1.) above (both with and without
;     Ricean corrections), though the folded analyses of (2.) and (3.) are both shown to the user.
;
;   DPC 20051026 original code
;       20051029 added stokes parameters calculations
;       20051226 fixed Q2U2 rotation matrix error (changed s -> -S) and redefined sinusoid PA (cos not sin)
;       20051229 found proper rotation amount is 45 deg, not 22.5
;       20060101 major rewrite to focus on Qi-Ui type analysis, as it returns most accurate uncertainties
;
file_in = dialog_pickfile(Title="Select Polarization Standard Star Summary File")
openr, lu_in, file_in, /GET_LUN
;
;   create structure for one star's aperture photometry data
;
star = {filename:'', hwp_angle:0.0, xcenter: fltarr(6), ycenter: fltarr(6), radii:fltarr(6), mag:dblarr(6), s_mag:dblarr(6), flux:dblarr(6), $
       s_flux:dblarr(6), nbad:intarr(6), s_noise:dblarr(4,6)}
;
;   and replicate this structure enough times for the polarimetry measurements (usually 32)
;
pol_data = REPLICATE(star,32)
;
;   read file
;
old_filename = ''
iap = 0
imaxap = 0
iframe = -1
imaxframe = 0
while (~EOF(lu_in)) do begin
    line = ''
    readf, lu_in, line
    parts = STRSPLIT(line,';',/EXTRACT)
    filelen = STRLEN(parts[0])
    rest = STRMID(line,filelen+1)       ; remove filename (only variable length field in input)
    ;
    ;  filename first
    ;
    if(parts[0] ne old_filename) then begin       ; new frame, reset counters
       ;
       old_filename = parts[0]
       iframe++
       imaxframe++
       iap=0
       ;
       pol_data[iframe].filename = parts[0]
       reads, rest, format='(f6.2)',HWP_angle
       pol_data[iframe].hwp_angle = HWP_angle
    endif
    ;
    ;  read remainder of information
    ;
    reads, rest, format='(7x,f6.1,1x,f6.1,1x,f6.1,1x,f9.5,1x,f8.5,1x,f13.2,1x,f9.2,1x,i3,4(1x,f9.1))',$
                  x,y,radius,mag,s_mag,flux,s_flux,nbad,s_noise_0,s_noise_1,$
                     s_noise_2,s_noise_3
    ;
    pol_data[iframe].xcenter[iap] = x
    pol_data[iframe].ycenter[iap] = y
    pol_data[iframe].radii[iap] = radius
    pol_data[iframe].mag[iap] = mag
    pol_data[iframe].s_mag[iap] = s_mag
    pol_data[iframe].flux[iap] = flux
    pol_data[iframe].s_flux[iap] = s_flux
    pol_data[iframe].nbad[iap] = nbad
    pol_data[iframe].s_noise[0,iap] = s_noise_0
    pol_data[iframe].s_noise[1,iap] = s_noise_1
    pol_data[iframe].s_noise[2,iap] = s_noise_2
    pol_data[iframe].s_noise[3,iap] = s_noise_3
    ;
    iap++
    imaxap = imaxap > iap
endwhile
;
FREE_LUN, lu_in
imaxframe--
pol_data = pol_data[0:imaxframe]
;
;   print summary
;
print, "Read in a total of ",imaxframe+1," frames and ",imaxap," apertures"
;
;for i = 0, imaxframe do begin
;    print, "frame ",i," HWP step ",pol_data[i].hwp_angle
;endfor
;
;   sort by HWP angle
;
sort_index = SORT(pol_data[*].hwp_angle)
pol_data = pol_data[sort_index]
;
;   scale HWP to angle
pol_data[*].hwp_angle = pol_data[*].hwp_angle * 0.675
;
;   plot all data for all apertures
;
x = pol_data[*].hwp_angle
xmax = MAX(x)
xmin = MIN(x)
;
dx = xmax - xmin
xmax = xmax + 0.05 * dx
xmin = xmin - 0.05 * dx
;
ymin = 1.0e6
for i = 0, imaxap-1 do begin
    ymin = ymin < MIN(pol_data[WHERE(pol_data[*].flux[i] ne 0.0)].flux[i])
endfor
ymax = MAX(pol_data[*].flux[*])
dy = ymax - ymin
ymax = ymax + 0.1 * dy
ymin = ymin - 0.1 * dy
;
;print, "Ymin, max = ",ymin, ymax
;
;   set up some colors
;
tvlct, [0, 255, 200, 0, 200, 0, 0, 255], [0, 0, 200, 0, 0, 255, 200, 255], [0, 0, 0, 255, 200, 0, 200, 255], 9
;
;   plot
;
device, decomposed=0
window, 11, Xsize=800, ysize=600, xpos=200, ypos=200
y = pol_data[WHERE(pol_data[*].flux[0] ne 0.0)].flux[0]
plot, x, y, XRANGE=[xmin,xmax], YRANGE=[ymin,ymax], Xstyle=1, thick=1, Ystyle=1, background=9, /NODATA, Xtitle="HWP Position Angle [deg.]",$
     Ytitle="Star Counts [ADUs]", Title="Standard Star Polarimetry"
;
for i = 0, imaxap-1 do begin
    y = pol_data[WHERE(pol_data[*].flux[i] ne 0.0)].flux[i]
    oplot, x, y, color=(i+10), thick=2, psym=-2
endfor
;
; slap up cursor and wait
;
;
cursor, xx, yy, /WAIT
;
wait, 1
;
;
wdelete, 11
;
;   ------------------------- Qi-Ui analysis starts -------------------
;
;   compute U, Q for each data pair
;
q = dblarr(6,2,4)
u = dblarr(6,2,4)
s_q = dblarr(6,2,4)
s_u = dblarr(6,2,4)
;
;   get the set of indices to the correct hwp angles
;
pa_ind = lonarr(32)
for iangle = 0, 31 do begin
    angle = iangle * 11.25
    ind = WHERE(abs(pol_data[*].hwp_angle - angle) lt 5.0,count)
    if(count gt 0) then pa_ind[iangle] = ind else pa_ind[iangle] = -99
    if(count gt 1) then print, "ERROR - multiple frames at HWP=",pol_data[ind].hwp_angle
endfor
;
;
for iap = 0, 5 do begin
    for ioffset = 0, 1 do begin    ; normal angle, and 45 deg set
       for i = 0, 3  do begin ; number of different samples
         ;
         ;    test for presence of correct data
         ;
         if(pa_ind[0+8*i+ioffset] ne -99 and pa_ind[4+8*i+ioffset] ne -99) then begin
          q[iap,ioffset,i] = -1.0*(-1.0^ioffset)*(pol_data[pa_ind[0+8*i+ioffset]].flux[iap] - pol_data[pa_ind[4+8*i+ioffset]].flux[iap])
          norm = abs(pol_data[pa_ind[0+8*i+ioffset]].flux[iap] + pol_data[pa_ind[4+8*i+ioffset]].flux[iap])
          q[iap,ioffset,i] = q[iap,ioffset,i] / norm
          s_q[iap,ioffset,i] = sqrt(pol_data[pa_ind[0+8*i+ioffset]].s_flux[iap]^2 + pol_data[pa_ind[4+8*i+ioffset]].s_flux[iap]^2)
          s_q[iap,ioffset,i] = (s_q[iap,ioffset,i]/norm)*sqrt(1.0+q[iap,ioffset,i]^2)
         endif else begin
          q[iap,ioffset,i] = 0.0
          s_q[iap,ioffset,i] = 1.0
         endelse
         if(pa_ind[2+8*i+ioffset] ne -99 and pa_ind[6+8*i+ioffset] ne -99) then begin
          u[iap,ioffset,i] = (-1.0^ioffset)*(pol_data[pa_ind[2+8*i+ioffset]].flux[iap] - pol_data[pa_ind[6+8*i+ioffset]].flux[iap])
          norm = abs(pol_data[pa_ind[2+8*i+ioffset]].flux[iap] + pol_data[pa_ind[6+8*i+ioffset]].flux[iap])
          u[iap,ioffset,i] = u[iap,ioffset,i] / norm
          s_u[iap,ioffset,i] = sqrt(pol_data[pa_ind[2+8*i+ioffset]].s_flux[iap]^2 + pol_data[pa_ind[6+8*i+ioffset]].s_flux[iap]^2)
          s_u[iap,ioffset,i] = (s_u[iap,ioffset,i]/norm)*sqrt(1.0d0 + u[iap,ioffset,i]^2)
         endif else begin
          u[iap,ioffset,i] = 0.0
          s_u[iap,ioffset,i] = 1.0
         endelse
       endfor
    endfor
endfor
;
;   rotate second set
;
;   coefficients for -45 deg rotation
;
c = 0.7071067812
s = -0.7071067812
;
for iap = 0, 5 do begin
    for i = 0, 3 do begin
        ;
        ;  test for missing data before rotating
        ;
        if(u[iap,1,i] ne 0.0 and q[iap,1,i] ne 0.0) then begin
            q_rot = -u[iap,1,i]*s + q[iap,1,i]*c
            u_rot = u[iap,1,i]*c + q[iap,1,i]*s
            s_q_rot = sqrt((c*s_q[iap,1,i])^2 + (s*s_u[iap,1,i])^2)
            s_u_rot = sqrt((c*s_u[iap,1,i])^2 + (s*s_q[iap,1,i])^2)
            ;
            q[iap,1,i] = q_rot
            u[iap,1,i] = u_rot
            s_q[iap,1,i] = s_q_rot
            s_u[iap,1,i] = s_u_rot
        endif else begin
            q[iap,1,i] = 0.0
            u[iap,1,i] = 0.0
            s_q[iap,1,i] = 1.
            s_u[iap,1,i] = 1.
        endelse
    endfor
endfor
;
;   print, Ui Qi values
;
for iap = 0, 5 do begin
    print, "for Aperture ",iap
    ;
    for i = 0, 3 do begin
       print,"Q(",2*i,") = ",q[iap,0,i],"+/-",s_q[iap,0,i],"  Q(",2*i+1,") = ",q[iap,1,i],"+/-",s_q[iap,1,i]
    endfor
    for i = 0, 3 do begin
       print,"U(",2*i,") = ",u[iap,0,i],"+/-",s_u[iap,0,i],"  U(",2*i+1,") = ",u[iap,1,i],"+/-",s_u[iap,1,i]
    endfor
endfor
;
;   compute mean U, Q from samples
;
q_low = 0.0
s_q_low = 1.0
u_low = 0.0
s_u_low = 1.0
;
for iap = 0, 5 do begin
    Q_vec = dblarr(8)
    U_vec = dblarr(8)
    sQ_vec = dblarr(8)
    sU_vec = dblarr(8)
    ;
    Q_vec[0:3] = REFORM(q[iap,0,*])
    Q_vec[4:7] = REFORM(q[iap,1,*])
    U_vec[0:3] = REFORM(u[iap,0,*])
    U_vec[4:7] = REFORM(u[iap,1,*])
    ;
    sQ_vec[0:3] = REFORM(s_q[iap,0,*])
    sQ_vec[4:7] = REFORM(s_q[iap,1,*])
    sU_vec[0:3] = REFORM(s_u[iap,0,*])
    sU_vec[4:7] = REFORM(s_u[iap,1,*])
    ;
    ;  means, uncertainties
    ;
    ind_q = WHERE(Q_vec ne 0.0,count_q)
    ind_u = WHERE(U_vec ne 0.0,count_u)
    ;
    if(count_q gt 0) then begin
        good_Q_vec = Q_vec(ind_q)
        zz_q = median_filtered_mean(good_Q_vec)
        ;
        ;   trap uncertainties to no less than propagated photometric...
        ;
        good_s_Q_vec = sQ_vec(ind_q)
        zz_sq = median_filtered_mean(good_s_Q_vec)
        q_err = zz_sq[0] / sqrt(count_q - 1.)
        zz_q[1] = q_err > zz_q[1]     ; larger of measured or predicted
    endif else begin
        zz_q = fltarr(2)
        zz_q[0] = 0.0
        zz_q[1] = 1.0
    endelse
    if(count_u gt 0) then begin
        good_U_vec = U_vec(ind_u)
        zz_u = median_filtered_mean(good_U_vec)
        ;
        good_s_U_vec = sU_vec(ind_u)
        zz_su = median_filtered_mean(good_s_U_vec)
        u_err = zz_su[0] / sqrt(count_u - 1.)
        zz_u[1] = u_err > zz_u[1]
    endif else begin
        zz_u = fltarr(2)
        zz_u[0] = 0.0
        zz_u[1] = 1.0
    endelse
    ;
    zz_q = zz_q * 100.0
    zz_u = zz_u * 100.0
    ;
    zz_q[1] = zz_q[1]/sqrt(count_q-1)
    zz_u[1] = zz_u[1]/sqrt(count_u-1)
    ;
    if(zz_q[1] lt s_q_low) then begin
        s_q_low = zz_q[1]
        q_low = zz_q[0]
    endif
    if(zz_u[1] lt s_u_low) then begin
        s_u_low = zz_u[1]
        u_low = zz_u[0]
    endif
    ;
    print, "Ap = ",iap,"  Mean Q = ",zz_q[0],"+/-",zz_q[1],"  U = ",zz_u[0],"+/-",zz_u[1]
    Puq = sqrt(zz_q[0]^2 + zz_u[0]^2)
    PAuq = (28.6478898d0 * atan(zz_u[0],zz_q[0])) MOD 180.
    if(PAuq lt 0.0) then PAuq = (PAuq + 360.0d0) MOD 180.
    ;
    ;  and their uncertainties
    ;
    sPuq = (1.0d0/Puq) * sqrt((zz_q[1]*zz_q[0])^2 + (zz_u[1]*zz_u[0])^2)
    sPAuq = 28.6478898d0 * (sPuq / Puq)
    print, format='("QiUi   Analysis: P= ",f9.4," +/- ",f9.4," %, Ang= ",f9.4," +/- ",f9.4," deg.")',Puq, sPuq,PAuq,sPAuq
    ;
endfor
;
;   and lowest uncertainty values
;
print, "Lowest Q, U uncertainty values: Q = ",q_low,"+/-",s_q_low,"  U = ",u_low,"+/-",s_u_low
Puq = sqrt(q_low^2 + u_low^2)
PAuq = (28.6478898d0 * atan(u_low,q_low)) MOD 180.
if(PAuq lt 0.0) then PAuq = (PAuq + 360.0d0) MOD 180.
sPuq = (1.0d0/Puq) * sqrt((s_q_low*q_low)^2 + (s_u_low*u_low)^2)
sPAuq = 28.6478898d0 * (sPuq / Puq)
print, format='("Min Unc. Analysis: P= ",f9.4," +/- ",f9.4," %, Ang= ",f9.4," +/- ",f9.4," deg.")',Puq, sPuq,PAuq,sPAuq
;
;   Ricean correction
;
corr = 1.0d0 - (sPuq/Puq)^2
if(corr gt 0.) then begin
    Pr = Puq * sqrt(1.0d0 - (sPuq/Puq)^2)
    sPAr = 28.6478898d0 * (sPuq / Pr)
endif else begin
    Pr = 0.0
    sPAr = 90.
endelse
print, "   Ricean corrected -- P = ",Pr,"+/-",sPuq,"%  Ang= ",PAuq,"+/-",sPAr," deg."
;
;-------------------------------- Qi-Ui Analysis End -------------------------------------------
;
;-------------------------------- Folded Flux Analyses Start ---------------------------------
;
;   plot the individual fluxes, coloring the different apertures differently
;
;   plot as symbol, with error bars
;
;   next, fold
;
pol_data[*].hwp_angle = (pol_data[*].hwp_angle*4) MOD 360.  ; make pol phase, mod 360
;
;   pick off the 359 angle stuff, wrap back around 0
;
ind_neg = WHERE(pol_data[*].hwp_angle gt 350,count)
if(count ne 0) then pol_data[ind_neg].hwp_angle = pol_data[ind_neg].hwp_angle - 360.
;
;   sort again
;
sort_index = SORT(pol_data[*].hwp_angle)
pol_data = pol_data[sort_index]
;for i = 0, imaxframe do begin
;    print, "frame ",i," pol angle = ",pol_data[i].hwp_angle
;endfor
;
;   plot folded data
;
xmin = -10.
xmax = 370.
x = pol_data[*].hwp_angle
;window, 12, Xsize=800, ysize=600, Xpos=0, ypos=0
y = pol_data[WHERE(pol_data[*].flux[0] ne 0.0)].flux[0]
;plot, x, y, XRANGE=[xmin,xmax], YRANGE=[ymin,ymax], Xstyle=1, thick=1, Ystyle=1, background=9, /NODATA, Xtitle="Polarization Position Angle [deg.]",$
;     Ytitle="Star Counts [ADUs]", Title="Folded"
;
for i = 0, imaxap-1 do begin
    y = pol_data[WHERE(pol_data[*].flux[i] ne 0.0)].flux[i]
;    oplot, x, y, color=(i+10), thick=1, psym=2
endfor
;
;----------------------
;
;   bin into 8 bins, computing mean, rms in each bin (*** with median_filtered_mean ??? ***)
;
y8 = fltarr(imaxap,8)
s_y8 = fltarr(imaxap,8)
x8 = float(indgen(8)*45.)
;
zz = dblarr(2)
for i = 0, 7 do begin
   angle = x8(i)
   ind = WHERE(abs(pol_data[*].hwp_angle - angle) lt 10.,count) ; select out proper data
   if(count ne 0) then begin
        ;
        ;    select data for each aperture
        ;
        for j = 0, imaxap-1 do begin
               vector = pol_data[ind].flux[j]
               noise_vector = pol_data[ind].s_flux[j]
               ;
               ;  compute cleaned mean, stdev
               ;
               if(count gt 1) then zz = median_filtered_mean(vector) else begin
                 zz[0] = vector
                 zz[1] = noise_vector
               endelse
               y8[j,i] = zz[0]
               s_y8[j,i] = zz[1]
               ;
               rms_save = zz[1]
               ;
               ;print, "for aperture ",j," angle ",angle," initial mean, rms = ",zz[0],zz[1]
               ;
               ;    for the values within 2.5 sigma, average them with weighting by their uncertainties
               ;
               if(count gt 1) then begin
                    w_ind = WHERE(abs(vector[*]-zz[0])/zz[1] lt 2.5,new_count)
                    ;print, "new_count = ",new_count
                    if(new_count gt 0) then begin
                        zmean = 0.0d0
                        zweight = 0.0d0
                        znum = 0L
                        for iave=0,new_count-1 do begin
                            wt = 1.0d0/(noise_vector[w_ind[iave]]^2)
                            zmean = zmean + vector[w_ind[iave]]*wt
                            zweight = zweight + wt
                            znum++
                        endfor
                        zmean = zmean / zweight
                        zrms = 1.0d0 / sqrt(zweight)
                        ;
                        y8[j,i] = zmean
                        s_y8[j,i] = zrms     ;- this probably understimates true error, including sky noise
                        ;
                        ;print, "   mean, old_uncertainty = ",zmean, zrms
                        ;
                        ;   try rms of values, reduced by 1/sqrt(N)
                        ;
                        s_y8[j,i] = rms_save/sqrt(new_count-1)
                        ;                       ;
                        ;print, "      new_uncertainty = ",s_y8[j,i]
                        ;
                    endif
                endif
               ;
        endfor
    endif else begin
        y8[j,i] = 0.0
        s_y8[j,i] = 10000000.0
    endelse
endfor
;
;   new plot of data
;
xmin = -10.
xmax = 370.
x = [x8,360.]
window, 13, Xsize=800, ysize=600
y = [REFORM(y8[0,*]),y8[0,0]]
plot, x, y, XRANGE=[xmin,xmax], YRANGE=[ymin,ymax], Xstyle=1, thick=1, Ystyle=1, background=9, /NODATA, Xtitle="Polarization Position Angle [deg.]",$
     Ytitle="Star Counts [ADUs]", Title="Folded, Averaged"
;
for j = 0, imaxap-1 do begin
    y = [REFORM(y8[j,*]),y8[j,0]]
    oplot, x, y, color=(j+10), thick=1, psym=2
endfor
;
;   print table of folded data
;
;print, "Angle                  Aperture Folded Counts / Uncertainties"
;print, "[deg]                         0            1               2                3                    4                      5"
;for i8 = 0, 7 do begin
;    print, x8[i8],y8[0,i8],y8[1,i8],y8[2,i8],y8[3,i8],y8[4,i8],y8[5,i8]
;    print, x8[i8],s_y8[0,i8],s_y8[1,i8],s_y8[2,i8],s_y8[3,i8],s_y8[4,i8],s_y8[5,i8]
;endfor
;
;   fit for each aperture
;
stokes_q = dblarr(imaxap)
stokes_u = dblarr(imaxap)
stokes_s_q = dblarr(imaxap)
stokes_s_u = dblarr(imaxap)
stokes_i_q_1 = 0.0d0
stokes_i_q_2 = 0.0d0
stokes_i_u_1 = 0.0d0
stokes_i_u_2 = 0.0d0
stokes_p = dblarr(imaxap)
stokes_s_p = dblarr(imaxap)
stokes_ang = dblarr(imaxap)
stokes_s_ang = dblarr(imaxap)
stokes_inten = dblarr(8)
stokes_s_inten = dblarr(8)
stokes_s_ND_q_1 = 0.0d0
stokes_s_ND_q_2 = 0.0d0
stokes_s_ND_u_1 = 0.0d0
stokes_s_ND_u_2 = 0.0d0
stokes_s_q_1 = 0.0d0
stokes_s_q_2 = 0.0d0
stokes_s_u_1 = 0.0d0
stokes_s_u_2 = 0.0d0
stokes_q_1 = 0.0d0
stokes_q_2 = 0.0d0
stokes_u_1 = 0.0d0
stokes_u_2 = 0.0d0
;
p = dblarr(imaxap)
sp = dblarr(imaxap)
pa = dblarr(imaxap)
spa = dblarr(imaxap)
;
;
;   coefficients for -45 deg rotation
;
c = 0.7071067812
s = -0.7071067812
;
for j = 0, imaxap-1 do begin
    ;---------------------------------------------------------------------------------------------------------
    ;  Stokes Parameters Analysis first:
    ;
    ;  load intensities
    ;
    stokes_inten = REFORM(y8[j,*])
    stokes_s_inten = REFORM(s_y8[j,*])
    ;
    ;  compute odd, even total counts
    ;
    stokes_i_q_1 = stokes_inten[0] + stokes_inten[4]
    stokes_i_q_2 = stokes_inten[1] + stokes_inten[5]
    stokes_i_u_1 = stokes_inten[2] + stokes_inten[6]
    stokes_i_u_2 = stokes_inten[3] + stokes_inten[7]
    ;
    ;print, "Total Intensities: q_1, q_2, u_1, u_2 = ",stokes_i_q_1, stokes_i_q_2, stokes_i_u_1, stokes_i_u_2
    ;
    ;  and stokes u, q (negate U to make angle sense correct on sky)
    ;
    stokes_q_1 = (stokes_inten[0] - stokes_inten[4])/stokes_i_q_1
    stokes_q_2 = (stokes_inten[1] - stokes_inten[5])/stokes_i_q_2
    stokes_u_1 = -(stokes_inten[2] - stokes_inten[6])/stokes_i_u_1
    stokes_u_2 = -(stokes_inten[3] - stokes_inten[7])/stokes_i_u_2
    ;
    ;print, "Stokes q1, q2, u1, u2 = ",stokes_q_1, stokes_q_2, stokes_u_1, stokes_u_2
    ;print, "Polarizations q1u1, q2u2 = ",sqrt(stokes_q_1^2 + stokes_u_1^2),sqrt(stokes_q_2^2 + stokes_u_2^2)
    ;
    ;  uncertainties in numerator or denominator
    ;
    stokes_s_ND_q_1 = sqrt(stokes_s_inten[0]^2 + stokes_s_inten[4]^2)
    stokes_s_ND_q_2 = sqrt(stokes_s_inten[1]^2 + stokes_s_inten[5]^2)
    stokes_s_ND_u_1 = sqrt(stokes_s_inten[2]^2 + stokes_s_inten[6]^2)
    stokes_s_ND_u_2 = sqrt(stokes_s_inten[3]^2 + stokes_s_inten[7]^2)
    ;
    ;   and uncertainties in unaveraged q, u
    ;
    stokes_s_q_1 = (stokes_s_ND_q_1 / stokes_i_q_1) * sqrt(1.0d0 + stokes_q_1^2)
    stokes_s_q_2 = (stokes_s_ND_q_2 / stokes_i_q_2) * sqrt(1.0d0 + stokes_q_2^2)
    stokes_s_u_1 = (stokes_s_ND_u_1 / stokes_i_u_1) * sqrt(1.0d0 + stokes_u_1^2)
    stokes_s_u_2 = (stokes_s_ND_u_2 / stokes_i_u_2) * sqrt(1.0d0 + stokes_u_2^2)
    ;
    ;print, "Stokes Uncertainties sq1, sq2, su1, su2 = ",stokes_s_q_1, stokes_s_q_2, stokes_s_u_1, stokes_s_u_2
    ;
    ;   before averaging, must rotate the "2" set by 45 degrees then average with "1" set
    ;
    stokes_u_2_r = stokes_u_2 * c + stokes_q_2 * s
    stokes_q_2_r = -stokes_u_2 * s + stokes_q_2 * c
    ;
    ;print, "Rotated polarization q2r.u2r = ",sqrt(stokes_u_2_r^2 + stokes_q_2_r^2)
    ;print, "Compare Q - 1(unrotated) = ",stokes_q_1," 2 (after rotation) = ",stokes_q_2_r
    ;print, "Compare U - 1(unrotated) = ",stokes_u_1," 2 (after rotation) = ",stokes_u_2_r
    ;
    ;   and rotate uncertainties too, though this will overestimate them -- should calc covariance
    ;
    stokes_s_q_2_r = sqrt((c*stokes_s_q_2)^2 + (s*stokes_s_u_2)^2)
    stokes_s_u_2_r = sqrt((c*stokes_s_u_2)^2 + (s*stokes_s_q_2)^2)
    ;
    ;print, "Compare sQ - 1(unrotated) = ",stokes_s_q_1," 2 (after rotation) = ",stokes_s_q_2_r
    ;print, "Compare sU - 1(unrotated) = ",stokes_s_u_1," 2 (after rotation) = ",stokes_s_u_2_r
    ;
    ;   weighted means
    ;
    stokes_q[j] = ((stokes_q_1/(stokes_s_q_1^2)) + (stokes_q_2_r/(stokes_s_q_2_r^2))) / $
               (1./(stokes_s_q_1^2)+1./(stokes_s_q_2_r^2))
    stokes_u[j] = ((stokes_u_1/(stokes_s_u_1^2)) + (stokes_u_2_r/(stokes_s_u_2_r^2))) / $
               (1./(stokes_s_u_1^2)+1./(stokes_s_u_2_r^2))
    ;
    ;  propagate to q, u uncertainties *****
    ;
    stokes_s_q[j] = sqrt(1./(1./(stokes_s_q_1^2)+1./(stokes_s_q_2_r^2)))
    stokes_s_u[j] = sqrt(1./(1./(stokes_s_u_1^2)+1./(stokes_s_u_2_r^2)))
    ;
    ;  now, calculate polarization, position angle
    ;
    stokes_p[j] = sqrt(stokes_q[j]^2 + stokes_u[j]^2)
    stokes_ang[j] = (28.6478898d0 * atan(stokes_u[j],stokes_q[j])) MOD 180.
    if(stokes_ang[j] lt 0.0) then stokes_ang[j] = (stokes_ang[j] + 360.0d0) MOD 180.
    ;
    ;  and their uncertainties
    ;
    stokes_s_p[j] = (1.0d0/stokes_p[j]) * sqrt((stokes_s_q[j]*stokes_q[j])^2 + (stokes_s_u[j]*stokes_u[j])^2)
    stokes_s_ang[j] = 28.6478898d0 * (stokes_s_p[j] / stokes_p[j])
    ;
    ;  scale q, u, p to percentages
    ;
    stokes_q[j] = 100.0 * stokes_q[j]
    stokes_s_q[j] = 100.0 * stokes_s_q[j]
    stokes_u[j] = 100.0 * stokes_u[j]
    stokes_s_u[j] = 100.0 * stokes_s_u[j]
    stokes_p[j] = 100.0 * stokes_p[j]
    stokes_s_p[j] = 100.0 * stokes_s_p[j]
    ;----------------------------------------------------------------------------------------------
    ;  Sine wave fitting Anaysis second:
    ;
    ;   first, subtract mean
    ;
    mean = total(y8[j,*]) / 8.
    flux = REFORM(y8[j,*]) - mean
    ;
    ;   Define a vector of weights.
    ;
    weights = 1.0/(REFORM(s_y8[j,*])*REFORM(s_y8[j,*]))
    ;
    ;  Provide an initial guess of the function's parameters.
    ;
    A = [0.02*mean,0.]
    ;
    ;  Compute the parameters.
    yfit = CURVEFIT(x8, flux, weights, A, SIGMA, FUNCTION_NAME='gfunct')
    ;
    ;   swap angle sense for sky PA true
    ;
    ;print, "Sine Wave fit (raw) - Amplitude = ",a[0],"+/-",sigma[0],' Phase = ',a[1],"+/-",sigma[1]
    ;
    model_a = fltarr(2)
    model_a[1] = a[1]
    a[1] = a[1] - 90.0
    ;
    ;   compute polarization
    ;
    p[j] = a[0] / mean
    sp[j] = sigma[0] / mean
    ;
    ;   if polarization negative, flip sign and shift 180 deg
    ;
    p[j] = 100.0 * p[j]
    sp[j] = 100.0 * sp[j]
    if(p[j] lt 0.) then begin
        a[0] = -a[0]
        p[j] = -p[j]
        a[1] = (a[1] + 180.0) MOD 360.
        model_a[1] = (model_a[1] + 180.0) MOD 360.
    endif
    ;
    a[1] = a[1] MOD 360.
    if(a[1] lt 0.) then a[1] = a[1] + 360.
    pa[j] = (a[1]/2.) mod 180.
    spa[j] = sigma[1]/2.
    ;
    ;-----------------------------------------------------------------------------
    ;
    ;   print results
    ;
    print, '--- for Aperture ',j
    ;
    print, format='("Stokes Analysis  : Q= ",f9.4," +/- ",f9.4,", U= ",f9.4," +/- ",f9.4," %")',stokes_q[j],stokes_s_q[j],stokes_u[j],stokes_s_u[j]
    print, format='("                                     P= ",f9.4," +/- ",f9.4," %, Ang= ",f9.4," +/- ",f9.4," deg.")',stokes_p[j],stokes_s_p[j],stokes_ang[j],stokes_s_ang[j]
    print, format='("Sinusoid Analysis: P= ",f9.4," +/- ",f9.4," %, Ang= ",f9.4," +/- ",f9.4," deg.")',p[j], sp[j],pa[j],spa[j]
    ;
    ;   compute model points
    ;
    xmodel = findgen(190)*2.0-2.0
    ymodel = mean + a[0] * sin((model_a[1]+xmodel)*0.01745)
    ;
    oplot, xmodel, ymodel, color=(j+10)
    ;
    plabel = string(FORMAT="('P[%]=',f5.2,'+/-',f4.2,' PA=',f6.1,'+/-',f4.1)",p[j],sp[j],pa[j],spa[j])
    xyouts, 20, mean-0.5*a[0], plabel
    ;
endfor
;
cursor, xx, yy, /wait
wait, 1
;wdelete, 12
wdelete, 13
;
print, "Save results to file? (Y/N)"
iyn = ''
read, iyn
if(STRUPCASE(iyn) eq 'Y') then begin
    file_out = dialog_pickfile(Title="Select Output Summary File",file=FILE_BASENAME(file_in,"phot.dat")+'pol.dat')
    openw, lu_out, file_out, /GET_LUN
    ;
    ;   old output data for folded fitting and sinusoidal fitting - removed 1/1/2006 DPC
    ;
    ;printf, lu_out, format='("File=",a," Aperture Radii=",9(f6.0,1x))',file_in,pol_data[0].radii[*]
    ;for j = 0, imaxap-1 do begin
    ;   printf, lu_out, format='(2(f6.1,1h;),12(f9.4,1h;),15(f11.2,1h;),f11.4)',pol_data[0].xcenter[j],pol_data[0].ycenter[j],$
    ;           stokes_q[j],stokes_s_q[j],stokes_u[j],stokes_s_u[j],$
    ;            stokes_p[j],stokes_s_p[j], stokes_ang[j],stokes_s_ang[j],$
    ;            p[j],sp[j],pa[j],spa[j],y8[j,0],s_y8[j,0],y8[j,1],s_y8[j,1],$
    ;             y8[j,2],s_y8[j,2],y8[j,3],s_y8[j,3],y8[j,4],s_y8[j,4],y8[j,5],s_y8[j,5],y8[j,6],s_y8[j,6],$
    ;             y8[j,7],s_y8[j,7]
    ;endfor
    ;
    ;   new output, based on Qi-Ui analysis and minimal uncertainties - added 1/1/2006 DPC
    ;
    printf, lu_out, format='(2(f6.1,1h;),9(f9.4,1h;),f9.4)',pol_data[0].xcenter[0],pol_data[0].ycenter[0],$
               q_low,s_q_low,u_low,s_u_low,$
                Puq,sPuq, PAuq,sPAuq,$
                Pr,sPAr                  ; Ricean corrected P value and PA uncertainty value
    ;
    FREE_LUN, lu_out
endif
;
print, FORMAT="('done!')"
end



pro gfunct, pa, a, f, pder
    ss = sin(((pa+a[1]) MOD 360.0)*0.01745)
    f = a[0] * ss
    IF N_PARAMS() GE 4 THEN $
    pder = [[ss], [A[0] * 0.01745 * cos(((pa+a[1]) MOD 360.0)*0.01745)]]
END