; $Id: ooechomagnetic__f_polanalysisgeneral.pro,v 1.7 2006/02/17 23:09:25 kneller Exp $
;###############################################################################
;
;FILENAME:
;        ooEchoMagnetic::F_PolAnalysisGeneral.pro
;
;
;PURPOSE:
;   Calclulate magnetic components of NSE data.
;
;CATEGORY:
;   Preliminary analysis.
;
;METHODS IN FILE:
;   ooEchoMagnetic::F_PolAnalysisGeneral
;
;METHODS:
;
;
;
; AUTHOR:
; Larry Kneller
; NIST Center for Neutron Research
; 100 Bureau Drive, Gaithersburg, MD 20899
; United States
; kneller@nist.gov  301-975-8839
; Feb 2 2006
;
; LICENSE:
; The software in this file is written by an employee of
; National Institute of Standards and Technology
; as part of the DAVE software project.
;
; The DAVE software package is not subject to copyright protection
; and is in the public domain. It should be considered as an
; experimental neutron scattering data reduction, visualization, and
; analysis system. As such, the authors assume no responsibility
; whatsoever for its use, and make no guarantees, expressed or
; implied, about its quality, reliability, or any other
; characteristic. The use of certain trade names or commercial
; products does not imply any endorsement of a particular product,
; nor does it imply that the named product is necessarily the best
; product for the stated purpose. We would appreciate acknowledgment
; if the DAVE software is used or if the code in this file is
; included in another product.
;
;###############################################################################


Function ooEchoMagnetic::F_PolAnalysisGeneral,Verbose,$
                                                xvals=xvals,$
                                                plotdata=plotdata,$
                                                eplotdata=eplotdata,$
                                                image=image,$
                                                eimage=eimage,$
                                                theimageindex=theimageindex,$
                                                magneticImage=magneticImage,$
                                                emagneticImage=emagneticImage,$
                                                spinCoherentImage=spinCoherentImage,$
                                                espinCoherentImage=espinCoherentImage,$
                                                CoherentAndIsotopeIncoherentImage=CoherentAndIsotopeIncoherentImage,$
                                                eCoherentAndIsotopeIncoherentImage=eCoherentAndIsotopeIncoherentImage,$
                                                TotalImage=TotalImage,$
                                                eTotalImage=eTotalImage,$
                                                pixel=pixel,$
                                                magq1d=magq1d,$
                                                sq1dtotal=sq1dtotal,$
                                                sq1dmagnetic=sq1dmagnetic,$
                                                sq1dspinIncoherent=sq1dspinIncoherent,$
                                                sq1dCoherentAndIsotopeIncoherent=sq1dCoherentAndIsotopeIncoherent,$
                                                timeaverage=timeaverage



;
;NAME:
;        ooEchoMagnetic::F_PolAnalysisGeneral
;
;PURPOSE:
;               Calculate components of S(Q) for the magnetic systems.
;PARAMETERS:
;               Verbose             Flag for print the results.
;
;KEYWORDS:
;               xvals               Fourier times
;               plotdata            An array containing the 4 data sets to be plotted (averaged over entire image???)
;               image               Time averaged image of choice theImageIndex
;               eimage              error bars for the image
;               theimageindex       The choice of image
;                                       theimageindex = 4  ;magnetic
;                                       theimageindex = 6  ;spin incoherent
;                                       theimageindex = 8  ;coherent and isotope incoherent
;                                       theimageindex = 10 ;total (SHOULD EQUAL S(Q)!!!!)
;
;               pixel               Pixel choice for single pixel return values.
;                                   (For plot of Mag variable v. time at each pixel.)
;                                   Value of pixel should be (ysize*ycoord + xcoord), i.e. a 1-d index
;
;SPECIFIC IMAGE CHOICES:
;               magneticImage
;               emagneticImage
;               spinCoherentImage
;               espinCoherentImage
;               CoherentAndIsotopeIncoherentImage
;               eCoherentAndIsotopeIncoherentImage
;               TotalImage
;               eTotalImage




;CURRENTLY I PRODUCE THE IMAGE.
;IT SHOULD BE A SIMPLE EXPANSION TO GET THE 4 IMAGES, AND USE AN OUTSIDE METHOD TO DO THE
;TIME AVERAGING FOR EACH Q-ARC.
;
;BUT I NEED TO TELL S(Q) TO PLOT ALL 4 COMPONENTS.  CAN I DO THIS BY HAVING THE MAGNETIC CASE POP UP 4 TIMES
;IN THE EXISTING S(Q) PLOTTER????
;
;
;






	;	This calculates the different scattering cross sections from a xyz polarization
	;	measurement (waves PCounts and PCountsErr2) and puts the result in a
	;	wave "Sc". The background is neglected. It is assumed that the counts are
	;	arranged in the following order: P perpendicular to Q and horizontal (x), DOWN and UP;
	;	P vertical (z), DOWN and UP; P parallel to Q and horizontal (y), DOWN and UP.
	;	It is assumed that UP is counted with the flipper ON. General formulae are used,
	;	no assumptions about inelasticity or symmetric alignment are made.
	;	If Verbose is true, the results are printed in the history area.

	;	This version takes different beam polarisations in all 3 directions into account.
	;	Consequently, there is only one expression for the magnetic scattering, not two.
	;	Supply the beam polarisation as waves eff_BeamPol[ndet][3] and eff_BeamPolErr2[ndet][3] (xzy).
	;	gPolEff accounts for another global depolarisation due to the sample.
	;	The flipper efficiency is still a global.

    if n_elements(verbose) eq 0 then verbose = 0
    if n_elements(timeaverage) eq 0 then timeaverage = 0

;020706
;
;THERE ARE TESTING PROBLEMS WITH VERBOSE SET TO 1, SO ALWAYS SET TO 0 FOR NOW.

    verbose = 0
    if n_elements(pixel) eq 0 then pixel = -1
    if n_elements(magq1d) eq 0 then magq1d = 0 ; SET THIS TO 0,1,2,3,4 FOR EACH OF THE COMPONENTS????

;
;	Wave PAC = PCounts
PAC = dblarr(6,self.x_dim*self.y_dim,self.no_of_fourier_times)   ;SET VALUE IN LOOP BELOW.
dimSize = size(PAC)
;	ndet = dimSize(PAC,1)
ndet = dimSize[2];self.x_dim*self.y_dim
;	ntof = dimSize(PAC,2)
ntof = dimSize[3];self.no_of_fourier_times



;	Wave Err2 = PCountsErr2
Err2 = dblarr(6,ndet,ntof)
;	Wave eff_BPW = eff_BeamPol
eff_BPW = dblarr(ndet,3)
eff_BPW[*,0] = self.eff_BPW
eff_BPW[*,1] = self.eff_BPW
eff_BPW[*,2] = self.eff_BPW  ;THIS IS ONLY A SCALAR RIGHT NOW

;	Wave eff_BPE2 = eff_BeamPolErr2
eff_BPE2 = dblarr(ndet,3)
eff_BPE2 = 0.02*eff_BPW      ;I DON'T HAVE THESE VALUES RIGHT NOW.  Just assume 2%.

;020706
;THIS IS ERROR SQUARED
eff_BPE2 = (0.02*eff_BPW)^2

;	Variable Pol,PolErr,FlipEff,FlipEffErr
Pol = self.P
PolErr = 0.02*Pol
FlipEff = self.f
FlipEffErr = 0.02*self.f

;	Variable/g gPolEff,gPieff
gPolEff = 0.0
gPieff = 0.0
;	Variable/g gPolEffErr,gPieffErr
gPolEffErr = 0.0
gPieffErr = 0.0

;	String cmd
cmd = ''

;	Make/O/N=(14,ndet,ntof) Sc			;	values and squared errors
Sc = dblarr(14,ndet,ntof)
;	Make/N=(ndet,ntof) xx,yy,zz,aa,bb,cc
xx = dblarr(ndet,ntof)
yy = dblarr(ndet,ntof)
zz = dblarr(ndet,ntof)
aa = dblarr(ndet,ntof)
bb = dblarr(ndet,ntof)
cc = dblarr(ndet,ntof)
Sc = dblarr(14,ndet,ntof)
sxx = dblarr(ndet,ntof)
syy = dblarr(ndet,ntof)
szz = dblarr(ndet,ntof)

xx_tave = dblarr(ndet)
yy_tave = dblarr(ndet)
zz_tave = dblarr(ndet)
aa_tave = dblarr(ndet)
bb_tave = dblarr(ndet)
cc_tave = dblarr(ndet)
Sc_tave = dblarr(14,ndet)
sxx_tave = dblarr(ndet)
syy_tave = dblarr(ndet)
szz_tave = dblarr(ndet)




;
;	; 1. Provide the beam polarization and the flipper efficiency.
	if (gPolEff eq 0.0) then begin
		gPolEff = 1.0
		gPolEffErr = 0.0
	endif
	if (gPieff eq 0.0) then begin
		gPieff = 1
		gPieffErr = 0
	endif
	Pol = gPolEff
	PolErr = gPolEffErr		;	will be neglected
	FlipEff = gPieff
	FlipEffErr = gPieffErr

;	; 2. dummies
	Sc[0,*,*] = 0.0
	Sc[1,*,*] = 0.0
	Sc[2,*,*] = 0.0
	Sc[3,*,*] = 0.0

;	; 3. Preparation
        for k=0,ntof-1 do begin

                tauind = k

;02/07/06
;THE Err2 VALUES ARE THE ERROR BARS SQUARED!!!!!
;ADDING THE SQUARE OPERATIONS TODAY.

                pac[0,*,k]  = self->helpChirpIup_downXYZArray(self.point_to_down_x,tauind)

                Err2[0,*,k] = (self->helpChirpIup_downXYZArray_eb(self.point_to_down_x,tauind))^2

                pac[1,*,k]  = self->helpChirpIup_downXYZArray(self.point_to_up_x,tauind)
                Err2[1,*,k] = (self->helpChirpIup_downXYZArray_eb(self.point_to_up_x,tauind))^2

                pac[4,*,k]  = self->helpChirpIup_downXYZArray(self.point_to_down_y,tauind)
                Err2[4,*,k] = (self->helpChirpIup_downXYZArray_eb(self.point_to_down_y,tauind))^2

                pac[5,*,k]  = self->helpChirpIup_downXYZArray(self.point_to_up_y,tauind)
                Err2[5,*,k] = (self->helpChirpIup_downXYZArray_eb(self.point_to_up_y,tauind))^2

                pac[2,*,k]  = self->helpChirpIup_downXYZArray(self.point_to_down_z,tauind)
                Err2[2,*,k] = (self->helpChirpIup_downXYZArray_eb(self.point_to_down_z,tauind))^2

                pac[3,*,k]  = self->helpChirpIup_downXYZArray(self.point_to_up_z,tauind)
                Err2[3,*,k] = (self->helpChirpIup_downXYZArray_eb(self.point_to_up_z,tauind))^2


                xx[*,k]  = self.xyzupdownscale*(pac[1,*,k]-pac[0,*,k])/(self.p*self.eff_BPW*(self.f+1.0))
                sxx[*,k] = self.xyzupdownscale*sqrt(Err2[1,*,k]^2.0 + Err2[0,*,k]^2.0)/ $
                                                                (self.p*self.eff_BPW*(self.f+1.0))

                yy[*,k]  = self.xyzupdownscale*(pac[5,*,k]-pac[4,*,k])/(self.p*self.eff_BPW*(self.f+1.0))
                syy[*,k] = self.xyzupdownscale*sqrt(Err2[5,*,k]^2.0 + Err2[4,*,k]^2.0)/ $
                                                                (self.p*self.eff_BPW*(self.f+1.0))

                zz[*,k]  = self.xyzupdownscale*(pac[3,*,k]-pac[2,*,k])/(self.p*self.eff_BPW*(self.f+1.0))
                szz[*,k] = self.xyzupdownscale*sqrt(Err2[3,*,k]^2.0 + Err2[2,*,k]^2.0)/ $
                                                                (self.p*self.eff_BPW*(self.f+1.0))


                aa[*,k] = self.xyzupdownscale*(self.f*pac[0,*,k] + pac[1,*,k])/(self.f+1.0)
                bb[*,k] = self.xyzupdownscale*(self.f*pac[2,*,k] + pac[3,*,k])/(self.f+1.0)
                cc[*,k] = self.xyzupdownscale*(self.f*pac[4,*,k] + pac[5,*,k])/(self.f+1.0)



                ;SET UP CRUDE TIME AVERAGES
                xx_tave = xx_tave + xx[*,k]
                yy_tave = yy_tave + yy[*,k]
                zz_tave = zz_tave + zz[*,k]
                aa_tave = aa_tave + aa[*,k]
                bb_tave = bb_tave + bb[*,k]
                cc_tave = cc_tave + cc[*,k]



        endfor;k



;THIS WILL BE HANDLED OUTSIDE, BUT IT MAY BE USEFUL TO HANDLE HERE AS WELL,
;SINCE ANTONIO REQUESTED A
        if timeaverage ne 0 then begin
            print,'f_polanal... test time average'
            for k = 0,ntof-1 do begin
                ;CALCULATE TIME AVERAGES
                xx[*,k] = xx_tave/ntof
                yy[*,k] = yy_tave/ntof
                zz[*,k] = zz_tave/ntof
                aa[*,k] = aa_tave/ntof
                bb[*,k] = bb_tave/ntof
                cc[*,k] = cc_tave/ntof
            endfor;k
        endif;timeaverage
;
;	; 4. Calculate the magnetic scattering
    Sc[4,*,*] = 4.0*zz-2.0*xx-2.0*yy


	Sc[5,*,*] = FlipEffErr*FlipEffErr*(-4*zz/(FlipEff+1)+2*xx/(FlipEff+1)+2*yy/(FlipEff+1))^2
	Sc[5,*,*] += eff_BPE2[0,0]*(2*xx/eff_BPW[0,0])^2
	Sc[5,*,*] += eff_BPE2[0,1]*(4*zz/eff_BPW[0,1])^2
	Sc[5,*,*] += eff_BPE2[0,2]*(2*yy/eff_BPW[0,2])^2
	Sc[5,*,*] += (Err2[3,*,*]+Err2[2,*,*])*(4/Pol/eff_BPW[0,1]/(FlipEff+1))^2		;	count err z
	Sc[5,*,*] += (Err2[1,*,*]+Err2[0,*,*])*(2/Pol/eff_BPW[0,0]/(FlipEff+1))^2		;	count err x
	Sc[5,*,*] += (Err2[5,*,*]+Err2[4,*,*])*(2/Pol/eff_BPW[0,2]/(FlipEff+1))^2		;	count err y





;
;	; 5. Calculate spin incoherent scattering

	Sc[6,*,*] = -9*zz+3*xx+3*yy+aa+bb+cc
	Sc[7,*,*] = FlipEffErr*FlipEffErr* $
	            ((9*zz-3*xx-3*yy-(Pol*eff_BPW[0,0]*xx+Pol*eff_BPW[0,1]*zz+Pol*eff_BPW[0,2]*yy))/(FlipEff+1))^2
	Sc[7,*,*] += eff_BPE2[0,0]*(3*xx/eff_BPW[0,0])^2
	Sc[7,*,*] += eff_BPE2[0,1]*(9*zz/eff_BPW[0,1])^2
	Sc[7,*,*] += eff_BPE2[0,2]*(3*yy/eff_BPW[0,2])^2
	Sc[7,*,*] += (Err2[1,*,*]*(Pol*eff_BPW[0,0]+3)^2+ $
	              Err2[0,*,*]*(FlipEff*Pol*eff_BPW[0,0]-3)^2)*(1/Pol/eff_BPW[0,0]/(FlipEff+1))^2		;	count err x
	Sc[7,*,*] += (Err2[5,*,*]*(Pol*eff_BPW[0,2]+3)^2+ $
	              Err2[4,*,*]*(FlipEff*Pol*eff_BPW[0,2]-3)^2)*(1/Pol/eff_BPW[0,2]/(FlipEff+1))^2		;	count err y
	Sc[7,*,*] += (Err2[3,*,*]*(Pol*eff_BPW[0,1]-9)^2+ $
	              Err2[2,*,*]*(FlipEff*Pol*eff_BPW[0,1]+9)^2)*(1/Pol/eff_BPW[0,1]/(FlipEff+1))^2		;	count err z
	Sc[6,*,*] /= 2
	Sc[7,*,*] /= 4




;	; 6. Calculate isotope incoherent and coherent scattering

	Sc[8,*,*] = 3*zz +3*xx +3*yy +aa +bb +cc
	Sc[9,*,*] = FlipEffErr*FlipEffErr*((-3*zz -3*xx -3*yy -(Pol*eff_BPW[0,0]*xx +Pol*eff_BPW[0,1]*zz +Pol*eff_BPW[0,2]*yy ))/(FlipEff+1))^2
	Sc[9,*,*] += eff_BPE2[0,0]*(3*xx /eff_BPW[0,0])^2
	Sc[9,*,*] += eff_BPE2[0,1]*(3*zz /eff_BPW[0,1])^2
	Sc[9,*,*] += eff_BPE2[0,2]*(3*yy /eff_BPW[0,2])^2
	Sc[9,*,*] += (Err2[1,*,*]*(Pol*eff_BPW[0,0]+3)^2+Err2[0,*,*]*(FlipEff*Pol*eff_BPW[0,0]-3)^2)*(1/Pol/eff_BPW[0,0]/(FlipEff+1))^2		;	count err x
	Sc[9,*,*] += (Err2[5,*,*]*(Pol*eff_BPW[0,2]+3)^2+Err2[4,*,*]*(FlipEff*Pol*eff_BPW[0,2]-3)^2)*(1/Pol/eff_BPW[0,2]/(FlipEff+1))^2		;	count err y
	Sc[9,*,*] += (Err2[3,*,*]*(Pol*eff_BPW[0,1]+3)^2+Err2[2,*,*]*(FlipEff*Pol*eff_BPW[0,1]-3)^2)*(1/Pol/eff_BPW[0,1]/(FlipEff+1))^2		;	count err z
	Sc[8,*,*] /= 6
	Sc[9,*,*] /= 36


;	; 7. Calculate the sum of all contributions

	Sc[10,*,*] = aa +bb +cc
	Sc[11,*,*] = FlipEffErr*FlipEffErr*((Pol*eff_BPW[0,0]*xx +Pol*eff_BPW[0,1]*zz +Pol*eff_BPW[0,2]*yy )/(FlipEff+1))^2
	Sc[11,*,*] += (Err2[0,*,*]+Err2[1,*,*]+Err2[2,*,*]+Err2[3,*,*]+Err2[4,*,*]+Err2[5,*,*])/(FlipEff+1)^2
	Sc[10,*,*] /= 1.5
	Sc[11,*,*] /= 2.25




;
;	; 8. Calculate the angle epsilon

    ww = dblarr(6,ndet,ntof)

	ww[5,*,*] = abs((PAC[0,*,*]-PAC[2,*,*])/(PAC[4,*,*]-PAC[2,*,*]))
	ww[1,*,*] = atan(sqrt(abs(PAC[0,*,*]-PAC[2,*,*])),sqrt(abs(PAC[4,*,*]-PAC[2,*,*])))
	ww[3,*,*] = 1/4/ww[5,*,*]/((1+ww[5,*,*])^2)/((PAC[4,*,*]-PAC[2,*,*])^4)*(Err2[0,*,*]*(PAC[4,*,*]-PAC[2,*,*])^2+Err2[2,*,*]*(PAC[0,*,*]-PAC[4,*,*])^2+Err2[4,*,*]*(PAC[0,*,*]-PAC[2,*,*])^2)
	ww[5,*,*] = abs((PAC[3,*,*]-PAC[1,*,*])/(PAC[3,*,*]-PAC[5,*,*]))
	ww[2,*,*] = atan(sqrt(abs(PAC[3,*,*]-PAC[1,*,*])),sqrt(abs(PAC[3,*,*]-PAC[5,*,*])))
	ww[4,*,*] = 1/4/ww[5,*,*]/((1+ww[5,*,*])^2)/((PAC[5,*,*]-PAC[3,*,*])^4)*(Err2[1,*,*]*(PAC[5,*,*]-PAC[3,*,*])^2+Err2[3,*,*]*(PAC[1,*,*]-PAC[5,*,*])^2+Err2[5,*,*]*(PAC[1,*,*]-PAC[3,*,*])^2)
	ww[0,*,*] = 1/ww[3,*,*]+1/ww[4,*,*]
	Sc[12,*,*] = 180/!PI*(ww[1,*,*]/ww[3,*,*]+ww[2,*,*]/ww[4,*,*])/ww[0,*,*]
	Sc[13,*,*] = 180/!PI*180/!PI/ww[0,*,*]
;	Sc[13,*,*] = 180/!PI*sqrt(1/ww[0,*,*])




;
;	; 9. Results - sum over the whole detector

;DO AVERAGE IN THE
;   1) T DIRECTION
;   2) x-y DIRECTION


;1)
sz = size(Sc)



;CREATE AREA AVERAGED- AND TIME AVERAGED- ARRAYS
m011_Sc = dblarr(14,1,ntof)
mtimeAv_Sc = dblarr(14,ndet,1)



wh = where(*self.mask2d gt 0,count)


;AREA AVERAGING
if count gt 0 then begin
    for j=0,ntof-1 do begin
        count = 0
;120705
;CONVERT THE NEXT TO ARRAY OPS
;        for i=0,ndet-1 do begin
;            if (*self.mask2d)[i,j] gt 0 then begin
;                for k=0,(sz[1]/2)-1 do begin
;                    ;if j eq 0 then print,2*k,'_____',Sc[2*k,i,j],'+/-',Sc[2*k+1,i,j]
;                    m011_Sc[2*k,0,j] += Sc[2*k,i,j]/Sc[2*k+1,i,j]^2.0
;                    m011_Sc[2*k+1,0,j] += 1.0d/Sc[2*k+1,i,j]^2.0
;                endfor;k
;                count++
;            endif
;        endfor;i

        wh_mask2d = where((*self.mask2d)[*,j] gt 0,count)

        if count gt 0 then begin
            for k=0,(sz[1]/2)-1 do begin

;ORIGINAL PRIOR TO 02/07/06
;                m011_Sc[2*k,0,j] = total((Sc[2*k,wh_mask2d,j]/Sc[2*k+1,wh_mask2d,j]^2.0),2)
;                m011_Sc[2*k+1,0,j] = total((1.0d/Sc[2*k+1,wh_mask2d,j]^2.0),2)

;I BELIEVE THAT THE ODD INDEX VALUES ARE THE ERROR BARS SQUARED!!!!
                m011_Sc[2*k,0,j] = total((Sc[2*k,wh_mask2d,j]/Sc[2*k+1,wh_mask2d,j]),2)
                m011_Sc[2*k+1,0,j] = total((1.0d/Sc[2*k+1,wh_mask2d,j]),2)



            endfor;k
        endif


        for k=0,(sz[1]/2)-1 do begin
;ORIGINAL PRIOR TO 02/07/06
;            m011_Sc[2*k,0,j] = m011_Sc[2*k,0,j]/m011_Sc[2*k+1,0,j];(m011_Sc[2*k+1,0,j] eq 0 ? double('NaN')  : m011_Sc[2*k,0,j]/m011_Sc[2*k+1,0,j])
;            m011_Sc[2*k+1,0,j] = 1.0d/sqrt(m011_Sc[2*k+1,0,j])

;THIS SHOULD BE UNCHANGED
            m011_Sc[2*k,0,j] = m011_Sc[2*k,0,j]/m011_Sc[2*k+1,0,j];(m011_Sc[2*k+1,0,j] eq 0 ? double('NaN')  : m011_Sc[2*k,0,j]/m011_Sc[2*k+1,0,j])
            m011_Sc[2*k+1,0,j] = 1.0d/sqrt(m011_Sc[2*k+1,0,j])
;            print,m011_Sc[2*k,0,j]
;            print,m011_Sc[2*k+1,0,j]
        endfor;k
    endfor;j
endif

;
;             IQT[kq,j] = IQT[kq,j] + (*self.S)[i,j]/((*self.sS)[i,j])^2
;            sIQT[kq,j] = sIQT[kq,j] + $
;                            1.0d/((*self.sS)[i,j])^2
;
;
;        IQT = IQT/sIQT
;        sIQT = 1.0d/sqrt(sIQT)


;TIME AVERAGING
;120705
;CONVERT TO ARRAY OPS IF POSSIBLE
if count gt 0 then begin
    for j=0,ndet-1 do begin
        count = 0
        for i=0,ntof-2 do begin
            if (*self.mask2d)[j,i] gt 0 then begin
                for k=0,sz[1]-1 do begin
                    mtimeAv_Sc[k,j,0] += Sc[k,j,i]
                endfor;k
                count++
            endif
        endfor;i
        for k=0,sz[1]-1 do begin
            mtimeAv_Sc[k,j,0] = (count eq 0 ? 0.0  : mtimeAv_Sc[k,j,0]/double(count))

        endfor;k
    endfor;j
endif




device,get_decomposed = dc
	if (Verbose ne 0) then begin

                    device,decomposed = 1
                    black = 0L & red = 255L & green = 256L*red & blue = 256L*green & white = red+green+blue

                    xvals = *self.fouriertime

            ;        plotdata = reform(m011_Sc[4,0,*],sz[3])
                    m011_Sc_all = [m011_Sc[4,0,*],m011_Sc[6,0,*],m011_Sc[8,0,*],m011_Sc[10,0,*]]

            ;IS THE NEXT LINE REALLY THE ERROR BARS OR SOMETHING FROM WHICH THEY WOULD BE DERIVED?????
                    em011_Sc_all = sqrt([m011_Sc[4+1,0,*],m011_Sc[6+1,0,*],m011_Sc[8+1,0,*],m011_Sc[10+1,0,*]])


                ;020606 SET UP INDIVIDUAL PIXEL SELECTION
                    ;HERE
                    if pixel ne -1 then begin
                    ;DO I HAVE THE CURRENTLY SELECTED PIXEL?
                        Sc_pixel = [Sc[4,pixel,*],Sc[6,pixel,*],Sc[8,pixel,*],Sc[10,pixel,*]]
                        eSc_pixel = sqrt([Sc[4+1,pixel,*],Sc[6+1,pixel,*],Sc[8+1,pixel,*],Sc[10+1,pixel,*]])

                        plotdata = reform(Sc_pixel,4*sz[3])
                        eplotdata = reform(eSc_pixel,4*sz[3])
                    endif else begin
                        plotdata = reform(m011_Sc_all,4*sz[3])
                        eplotdata = reform(em011_Sc_all,4*sz[3])
                    endelse


                    window,0,xsize=300,ysize=300,xpos=0,ypos=0,title='Component plot'
                    plot,*self.fouriertime,reform(m011_Sc[4,0,*],sz[3]),background = white,color=black,$;xrange = [0,30],$
                                                $;yrange=[min(m011_Sc[4,0,*]),max(m011_Sc[4,0,*])],$;;;yrange = [-700,4700],$[min(m011_Sc),max(m011_Sc)],$
                                                yrange=[min(m011_Sc_all),max(m011_Sc_all)],$
                                                title='Components: Area averaged',psym=-2,/xlog,$
                                                xtitle = 'Fourier Time'
                    oplot,*self.fouriertime,reform(m011_Sc[6,0,*],sz[3]),color=red,psym=-4
                    oplot,*self.fouriertime,reform(m011_Sc[8,0,*],sz[3]),color=green,psym=-5
                    oplot,*self.fouriertime,reform(m011_Sc[10,0,*],sz[3]),color=blue,psym=-6

                    xo=0.55 & yo = 0.5
                    nse_legend,0,xo,yo,'magnetic',color=black,psym=-2
                    nse_legend,0,xo,yo-0.04,'spin incoh',color=red,psym=-4
                    nse_legend,0,xo,yo-0.08,'coh & isotope incoh',color=green,psym=-5
                    nse_legend,0,xo,yo-0.12,'total',color=blue,psym=-6

                    device,decomposed=0
                    window,1,xsize=300,ysize=300,xpos = 0, ypos=300, title = 'Image'
                    ;print,max(mtimeAv_Sc[4,*,0])


                    tvlct,rnow,gnow,bnow,/get

                    ;EXPERIMENTING WITH SETTING BACKGROUND WHITE, FOREGROUND BLACK
                    rnownew = rnow
                    gnownew = gnow
                    bnownew = bnow

                    rnownew[0] = 255L
                    gnownew[0] = 255L
                    bnownew[0] = 255L

                    ;SET THE MAXIMUM COLOR TO BLACK
                    rnownew[255] = 0L
                    gnownew[255] = 0L
                    bnownew[255] = 0L
                    impos = [0.24,0.20,0.9,0.9]
                    cbpos = [0.20,0.13,0.9,0.16]

                    tvlct,rnownew,gnownew,bnownew

                    if n_elements(theimageindex) eq 0 then theimageindex = 0
                    theimageindex = 4+2*theimageindex
                    ;theimageindex = 4  ;magnetic
                    ;theimageindex = 6  ;spin incoherent
                    ;theimageindex = 8  ;coherent and isotope incoherent
                    ;theimageindex = 10 ;total
                    image = reform(mtimeAv_Sc[theimageindex,*,0],sqrt(ndet),sqrt(ndet))
                    eimage = reform(mtimeAv_Sc[theimageindex+1,*,0],sqrt(ndet),sqrt(ndet))

                    theim = bytscl(image,top=253)+1B

                    plotimage,theim,yrange = [8,0],$
                                        title='magnetic: time averaged',$
                                        position=impos
                                        ;bottom = white,$
                                        ;top = black

                    device,decomposed=1
                        for i=0,ndet-1 do begin
                            xpix = i mod self.x_dim
                            ypix = i/self.x_dim
                            if theIm[i] eq 0B then begin
                                x_xout = [xpix,xpix+1,xpix+1,xpix,xpix]
                                y_xout = [ypix,ypix,ypix+1,ypix+1,ypix]
                                polyfill,x_xout,y_xout,color=green;xcolor2

                            endif
                        endfor;i
                    device,decomposed=0


                    COLORBAR,   NCOLORS=256, $
                                $;POSITION=[0.10, 0.15, 0.15, 0.90],$
                                bottom=white,$
                                top =black,$
                                range = [min(mtimeAv_Sc[theimageindex,*,0]),max(mtimeAv_Sc[theimageindex,*,0])],$
                                /vertical,$
                                position=cbpos

                    ;RESET THE COLORS
                    tvlct,rnow,gnow,bnow
                    device,decomposed=1
                    device,decomposed = dc

                        print, format='(a,f7.3,a,f6.3,a)','GENERAL Pol. Analysis gives magnetic:',$
                            		                      double(m011_Sc[4,0,0]),' ± ',sqrt(double(m011_Sc[5,0,0])),$
                                                          ' Counts/sec (4*gPol/gPolEff)'


                        print, format='(a,f7.3,a,f6.3,a)','spin incoherent:',$
                                                            double(m011_Sc[6,0,0]),' ± ',sqrt(double(m011_Sc[7,0,0])),$
                                                            'Counts/sec'

                		print, format='(a,f7.3,a,f6.3,a)','coherent and isotope incoherent: ',$
                		                                    double(m011_Sc[8,0,0]),' ± ',sqrt(double(m011_Sc[9,0,0])),$
                		                                    ' Counts/sec'

                        print, format='(a,f7.3,a,f6.3,a)','sum of all contributions: ',$
                                                          double(m011_Sc[10,0,0]),' ± ',$
                                                          sqrt(double(m011_Sc[11,0,0])),' Counts/sec (2*gAver)'

                        print, format='(a,f7.3,a,f6.3,a)','angle epsilon: ',$
                                                          double(m011_Sc[12,0,0]),' ± ',$
                                                          sqrt(double(m011_Sc[13,0,0])),' (deg)'

    endif else begin
                    ;NOT VERBOSE

                    if n_elements(theimageindex) eq 0 then theimageindex = 0
                    theimageindex = 4+2*theimageindex
                    ;theimageindex = 4  ;magnetic
                    ;theimageindex = 6  ;spin incoherent
                    ;theimageindex = 8  ;coherent and isotope incoherent
                    ;theimageindex = 10 ;total
                    if arg_present(image) then $
                        image = reform(mtimeAv_Sc[theimageindex,*,0],sqrt(ndet),sqrt(ndet))
                    if arg_present(eimage) then $
                        eimage = reform(mtimeAv_Sc[theimageindex+1,*,0],sqrt(ndet),sqrt(ndet))


                ;020606 SET UP INDIVIDUAL PIXEL SELECTION
                    ;HERE

                    if arg_present(plotdata) gt 0 or arg_present(eplotdata) gt 0 then begin
                        if pixel ne -1 then begin
                        ;DO I HAVE THE CURRENTLY SELECTED PIXEL?
                            ;print,pixel
                            Sc_pixel = [Sc[4,pixel,*],Sc[6,pixel,*],Sc[8,pixel,*],Sc[10,pixel,*]]
                            eSc_pixel = sqrt([Sc[4+1,pixel,*],Sc[6+1,pixel,*],Sc[8+1,pixel,*],Sc[10+1,pixel,*]])
                            plotdata = reform(Sc_pixel,4*sz[3])
                            eplotdata = reform(eSc_pixel,4*sz[3])
                        endif else begin
                                ;plotdata = reform(m011_Sc[4,0,*],sz[3])
                                m011_Sc_all = [reform(m011_Sc[4,0,*],sz[3]),$
                                                reform(m011_Sc[6,0,*],sz[3]),$
                                                reform(m011_Sc[8,0,*],sz[3]),$
                                                reform(m011_Sc[10,0,*],sz[3])]


                        ;IS THE NEXT LINE REALLY THE ERROR BARS OR SOMETHING FROM WHICH THEY WOULD BE DERIVED?????
                                em011_Sc_all = sqrt([m011_Sc[4+1,0,*],m011_Sc[6+1,0,*],m011_Sc[8+1,0,*],m011_Sc[10+1,0,*]])
                            plotdata = reform(m011_Sc_all,4*sz[3])
                            eplotdata = reform(em011_Sc_all,4*sz[3])
                        endelse
                        ;help,m011_Sc_all
                        ;help,plotdata
                    endif;arg_present(plotdata)
    endelse
;#############################################

;TUESDAY!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


;HAVING THESE IMAGES AND THE ABILITY TO GET THE Q ARCS MATRIX AND Q VALUES MEANS THAT I CAN
;EASILY MAKE EACH COMPONENT OF THE S(Q) AVAILABLE.


;#############################################

    if arg_present(magneticImage) gt 0 then begin
            MagneticImage = reform(mtimeAv_Sc[4,*,0],sqrt(ndet),sqrt(ndet))
    endif
    if arg_present(emagneticImage) gt 0 then begin
            eMagneticImage = sqrt(reform(mtimeAv_Sc[4+1,*,0],sqrt(ndet),sqrt(ndet)))
    endif
    if arg_present(spinCoherentImage) gt 0 then begin
            spinCoherentImage = reform(mtimeAv_Sc[6,*,0],sqrt(ndet),sqrt(ndet))
    endif
    if arg_present(espinCoherentImage) gt 0 then begin
            espinCoherentImage = sqrt(reform(mtimeAv_Sc[6+1,*,0],sqrt(ndet),sqrt(ndet)))
    endif
    if arg_present(CoherentAndIsotopeIncoherentImage) gt 0 then begin
            CoherentAndIsotopeIncoherentImage = reform(mtimeAv_Sc[8,*,0],sqrt(ndet),sqrt(ndet))
    endif
    if arg_present(eCoherentAndIsotopeIncoherentImage) gt 0 then begin
            eCoherentAndIsotopeIncoherentImage = sqrt(reform(mtimeAv_Sc[8+1,*,0],sqrt(ndet),sqrt(ndet)))
    endif

if arg_present(SQ1DTotal) gt 0 or $
    arg_present(TotalImage) gt 0 or $
     arg_present(eTotalImage) gt 0 then begin

            ;mtimeAv_Sc SEEMS TO CORRECTLY ACCOUNT FOR selfmask2d ABOVE
            ;FOR THE DATA, BUT NOT FOR THE ERROR BARS!!!!!!!!!!!!!!!!!!


            ;;;;;
            ;;;;;;;JUST GET Q INDEX MATRIX, THE Q VALUES AND THE MASK CREATE THE TIME AVERAGED VALUES
            ;;;;;;
            ;;;;;        self->calculateArcsQMatrix,narcs=narcs,arcmat=arcmat
            ;;;;;;    ;THIS PROCEDURE INCLUDES self->calculateQVals WHICH CALCULATES THE VALUES OF Q
            ;;;;;;    ;AND STORES THEM IN self.Qvals
            ;;;;;    ;THE Q-VALUES OF INTEREST HERE ARE THE ONES IN self.qarcs

                mask = *(self->getproperty(tag='mask2d'))
                sz = size(Sc)
                whrmask = where(mask le 0,count)
                if count ne 0 then begin
                    mask[whrmask] = 0
                endif

                ;help,mask

                etot = reform(Sc[10,*,*],sz[2],sz[3])  ;total(e[*,0:self->term()-1,*],2)
                setot = reform(sqrt(Sc[11,*,*]),sz[2],sz[3])

                self->calculateArcsQMatrix,narcs=narcs,arcmat=mat

                if n_elements(narcs) eq 0 then narcs = self.narcs

                theQvals = *(self->getproperty(tag='Qvals'))

                s = dblarr(narcs)
                sS = dblarr(narcs)
                q = dblarr(narcs)
                qcount = intarr(narcs)
                sz = size(etot)
                for i=0,sz[1]-1 do begin
                    for j=0,sz[2]-1 do begin
                        for k=0,narcs-1 do begin
                            if mat[i] eq k and mask[i,j] ne 0 then begin
                                if finite (etot[i,j]) gt 0  and finite(1.0d/setot[i,j]) gt 0 then begin
                                    s[k] += etot[i,j]/setot[i,j]^2
                                    ss[k] += 1.0d/setot[i,j]^2

                                    q[k] += theQVals[i,j]
                                    qcount[k] += 1
                                endif;finite

                            endif

                        endfor;k
                    endfor;j
                endfor;i
                ;print,'qcount=',qcount
                qmin = min(theqvals)
                qmax = max(theqvals)

                for i=0,narcs-1 do begin
                    if qcount[i] ne 0 then begin
                        s[i] = s[i]/ss[i]
                        ss[i] = 1.0d/sqrt(ss[i])
                        q[i] = q[i]/double(qcount[i])
                    endif else begin
                        s[i] = 0.0
                        ss[i] = 0.0
                        if narcs ne 1 then begin
                            qstep = (qmax-qmin)/double(narcs-1)
                            q[i] = qmin + (double(i)+0.5)*qstep
                        endif else begin
                            q[i] = self->getproperty(tag='q')
                        endelse
                    endelse
                endfor;i
;print,self.preset
                s = s/double(self.preset)
;print,s

                ss = ss/double(self.preset)
                ;print,q
                ;print,s
                ;print,ss

                SQ1DTotal = {s:ptr_new(s),ss:ptr_new(ss),q:ptr_new(q)}
            ;    window,0
            ;    wset,0
            ;    plot,[q],[s],psym=-2
            ;    errplot,[q],[s-ss],[s+ss]

endif;SQ1DTotal


if arg_present(sq1dmagnetic) gt 0 or $
    arg_present(TotalImage) gt 0 or $
     arg_present(eTotalImage) gt 0 then begin

            ;mtimeAv_Sc SEEMS TO CORRECTLY ACCOUNT FOR selfmask2d ABOVE
            ;FOR THE DATA, BUT NOT FOR THE ERROR BARS!!!!!!!!!!!!!!!!!!


            ;;;;;
            ;;;;;;;JUST GET Q INDEX MATRIX, THE Q VALUES AND THE MASK CREATE THE TIME AVERAGED VALUES
            ;;;;;;
            ;;;;;        self->calculateArcsQMatrix,narcs=narcs,arcmat=arcmat
            ;;;;;;    ;THIS PROCEDURE INCLUDES self->calculateQVals WHICH CALCULATES THE VALUES OF Q
            ;;;;;;    ;AND STORES THEM IN self.Qvals
            ;;;;;    ;THE Q-VALUES OF INTEREST HERE ARE THE ONES IN self.qarcs

                mask = *(self->getproperty(tag='mask2d'))
                sz = size(Sc)
                whrmask = where(mask le 0,count)
                if count ne 0 then begin
                    mask[whrmask] = 0
                endif

                ;help,mask

                etot = reform(Sc[4,*,*],sz[2],sz[3])  ;total(e[*,0:self->term()-1,*],2)
                setot = reform(sqrt(Sc[5,*,*]),sz[2],sz[3])



                self->calculateArcsQMatrix,narcs=narcs,arcmat=mat

                if n_elements(narcs) eq 0 then narcs = self.narcs

                theQvals = *(self->getproperty(tag='Qvals'))

                s = dblarr(narcs)
                sS = dblarr(narcs)
                q = dblarr(narcs)
                qcount = intarr(narcs)

                sz = size(etot)
                for i=0,sz[1]-1 do begin
                    for j=0,sz[2]-1 do begin
                        for k=0,narcs-1 do begin
                            if mat[i] eq k and mask[i,j] ne 0 then begin
                                if finite (etot[i,j]) gt 0  and finite(1.0d/setot[i,j]) gt 0 then begin
                                    s[k] += etot[i,j]/setot[i,j]^2
                                    ss[k] += 1.0d/setot[i,j]^2

                                    q[k] += theQVals[i,j]
                                    qcount[k] += 1
                                endif;finite
                            endif
                        endfor;k
                    endfor;j
                endfor;i

                qmin = min(theqvals)
                qmax = max(theqvals)

                for i=0,narcs-1 do begin
                    if qcount[i] ne 0 then begin
                        s[i] = s[i]/ss[i]
                        ss[i] = 1.0d/sqrt(ss[i])
                        q[i] = q[i]/double(qcount[i])
                    endif else begin
                        s[i] = 0.0
                        ss[i] = 0.0
                        if narcs ne 1 then begin
                            qstep = (qmax-qmin)/double(narcs-1)
                            q[i] = qmin + (double(i)+0.5)*qstep
                        endif else begin
                            q[i] = self->getproperty(tag='q')
                        endelse
                    endelse
                endfor;i

                s = s/double(self.preset)
                ss = ss/double(self.preset)
                ;print,q
                ;print,s
                ;print,ss

                sq1dmagnetic = {s:ptr_new(s),ss:ptr_new(ss),q:ptr_new(q)}
            ;    window,0
            ;    wset,0
            ;    plot,[q],[s],psym=-2
            ;    errplot,[q],[s-ss],[s+ss]

endif;sq1dmagnetic






if arg_present(sq1dspinIncoherent) gt 0 or $
    arg_present(TotalImage) gt 0 or $
     arg_present(eTotalImage) gt 0 then begin

            ;mtimeAv_Sc SEEMS TO CORRECTLY ACCOUNT FOR selfmask2d ABOVE
            ;FOR THE DATA, BUT NOT FOR THE ERROR BARS!!!!!!!!!!!!!!!!!!


            ;;;;;
            ;;;;;;;JUST GET Q INDEX MATRIX, THE Q VALUES AND THE MASK CREATE THE TIME AVERAGED VALUES
            ;;;;;;
            ;;;;;        self->calculateArcsQMatrix,narcs=narcs,arcmat=arcmat
            ;;;;;;    ;THIS PROCEDURE INCLUDES self->calculateQVals WHICH CALCULATES THE VALUES OF Q
            ;;;;;;    ;AND STORES THEM IN self.Qvals
            ;;;;;    ;THE Q-VALUES OF INTEREST HERE ARE THE ONES IN self.qarcs

                mask = *(self->getproperty(tag='mask2d'))
                sz = size(Sc)
                whrmask = where(mask le 0,count)
                if count ne 0 then begin
                    mask[whrmask] = 0
                endif

                ;help,mask

                etot = reform(Sc[6,*,*],sz[2],sz[3])  ;total(e[*,0:self->term()-1,*],2)
                setot = reform(sqrt(Sc[7,*,*]),sz[2],sz[3])



                self->calculateArcsQMatrix,narcs=narcs,arcmat=mat

                if n_elements(narcs) eq 0 then narcs = self.narcs

                theQvals = *(self->getproperty(tag='Qvals'))

                s = dblarr(narcs)
                sS = dblarr(narcs)
                q = dblarr(narcs)
                qcount = intarr(narcs)

                sz = size(etot)
                for i=0,sz[1]-1 do begin
                    for j=0,sz[2]-1 do begin
                        for k=0,narcs-1 do begin
                            if mat[i] eq k and mask[i,j] ne 0 then begin
                                if finite (etot[i,j]) gt 0  and finite(1.0d/setot[i,j]) gt 0 then begin

                                    s[k] += etot[i,j]/setot[i,j]^2
                                    ss[k] += 1.0d/setot[i,j]^2

                                    q[k] += theQVals[i,j]
                                    qcount[k] += 1
                                endif;finite
                            endif
                        endfor;k
                    endfor;j
                endfor;i

                qmin = min(theqvals)
                qmax = max(theqvals)

                for i=0,narcs-1 do begin
                    if qcount[i] ne 0 then begin
                        s[i] = s[i]/ss[i]
                        ss[i] = 1.0d/sqrt(ss[i])
                        q[i] = q[i]/double(qcount[i])
                    endif else begin
                        s[i] = 0.0
                        ss[i] = 0.0
                        if narcs ne 1 then begin
                            qstep = (qmax-qmin)/double(narcs-1)
                            q[i] = qmin + (double(i)+0.5)*qstep
                        endif else begin
                            q[i] = self->getproperty(tag='q')
                        endelse
                    endelse
                endfor;i

                s = s/double(self.preset)
                ss = ss/double(self.preset)
                ;print,q
                ;print,s
                ;print,ss

                sq1dspinIncoherent = {s:ptr_new(s),ss:ptr_new(ss),q:ptr_new(q)}
            ;    window,0
            ;    wset,0
            ;    plot,[q],[s],psym=-2
            ;    errplot,[q],[s-ss],[s+ss]

endif;sq1dspinIncoherent


if arg_present(sq1dCoherentAndIsotopeIncoherent) gt 0 or $
    arg_present(TotalImage) gt 0 or $
     arg_present(eTotalImage) gt 0 then begin

            ;mtimeAv_Sc SEEMS TO CORRECTLY ACCOUNT FOR selfmask2d ABOVE
            ;FOR THE DATA, BUT NOT FOR THE ERROR BARS!!!!!!!!!!!!!!!!!!


            ;;;;;
            ;;;;;;;JUST GET Q INDEX MATRIX, THE Q VALUES AND THE MASK CREATE THE TIME AVERAGED VALUES
            ;;;;;;
            ;;;;;        self->calculateArcsQMatrix,narcs=narcs,arcmat=arcmat
            ;;;;;;    ;THIS PROCEDURE INCLUDES self->calculateQVals WHICH CALCULATES THE VALUES OF Q
            ;;;;;;    ;AND STORES THEM IN self.Qvals
            ;;;;;    ;THE Q-VALUES OF INTEREST HERE ARE THE ONES IN self.qarcs

                mask = *(self->getproperty(tag='mask2d'))
                sz = size(Sc)
                whrmask = where(mask le 0,count)
                if count ne 0 then begin
                    mask[whrmask] = 0
                endif

                ;help,mask

                etot = reform(Sc[8,*,*],sz[2],sz[3])  ;total(e[*,0:self->term()-1,*],2)
;                help,etot
;                print,'ndet=',ndet
;                print,'ntof=',ntof
                setot = reform(sqrt(Sc[9,*,*]),sz[2],sz[3])
;print,'coh_____________________________'
;help,'finite(etot)=';,finite(etot)
;help,'where(finite(etot))=';,where(finite(etot))
;help,'where(finite(setot))=';,where(finite(setot))
;help,max(etot)
;help,max(setot)
;help,min(etot)
;help,min(setot)
                self->calculateArcsQMatrix,narcs=narcs,arcmat=mat

                if n_elements(narcs) eq 0 then narcs = self.narcs

                theQvals = *(self->getproperty(tag='Qvals'))

                s = dblarr(narcs)
                sS = dblarr(narcs)
                q = dblarr(narcs)
                qcount = intarr(narcs)

                sz = size(etot)


;print,size(etot)
;print,size(setot)
;print,min(mask)
                for i=0,sz[1]-1 do begin
                    for j=0,sz[2]-1 do begin
                        for k=0,narcs-1 do begin
                            if mat[i] eq k and mask[i,j] ne 0 then begin
                                if finite (etot[i,j]) gt 0  and finite(1.0d/setot[i,j]) gt 0 then begin
                                    s[k] += etot[i,j]/setot[i,j]^2
                                    ss[k] += 1.0d/setot[i,j]^2
                                    ;if finite(s[k]) eq 0 then print,i,j,k
                                    q[k] += theQVals[i,j]
                                    qcount[k] += 1
                                endif;finite
                            endif
                        endfor;k
                    endfor;j
                endfor;i
;print,'qcount=',qcount
                qmin = min(theqvals)
                qmax = max(theqvals)

                for i=0,narcs-1 do begin
                    if qcount[i] ne 0 then begin
                        s[i] = s[i]/ss[i]
                        ss[i] = 1.0d/sqrt(ss[i])
                        q[i] = q[i]/double(qcount[i])
                    endif else begin
                        s[i] = 0.0
                        ss[i] = 0.0
                        if narcs ne 1 then begin
                            qstep = (qmax-qmin)/double(narcs-1)
                            q[i] = qmin + (double(i)+0.5)*qstep
                        endif else begin
                            q[i] = self->getproperty(tag='q')
                        endelse
                    endelse
                endfor;i

                s = s/double(self.preset)
                ss = ss/double(self.preset)
                ;print,q
                ;print,s
                ;print,ss

                sq1dCoherentAndIsotopeIncoherent = {s:ptr_new(s),ss:ptr_new(ss),q:ptr_new(q)}
            ;    window,0
            ;    wset,0
            ;    plot,[q],[s],psym=-2
            ;    errplot,[q],[s-ss],[s+ss]

endif;sq1dCoherentAndIsotopeIncoherent





    return,1
end;ooechomagnetic::f_polanalysisgeneral
