;+ ;Name: ; spinfit ;Purpose: ; performs a spinfit on B or E field data, results should be ; equivalent to FGS or EFS datatypes ;Calling Sqeuence: ; spinfit,arr_in_t,arr_in_data,arr_in_sunpulse_t,arr_in_sunpulse_data, $ ; A,B,C, avg_axis, median_axis,Sigma,Npoints,sun_data, $ ; min_points=min_points,alpha=alpha,beta=beta, $ ; plane_dim=plane_dim,axis_dim=axis_dim,phase_mask_starts=phase_mask_starts, $ ; phase_mask_ends=phase_mask_ends,sun2sensor=sun2sensor ;Input: ; arr_in_t = time array for the data ; arr_in_data = the data to be spin fit ; arr_in_sunpulse_t = time array for sunpulse data ; arr_in_sunpulse_data = sunpulse data ;Output: ; A,B,C = fit parameters for spinfit ; avg_axis = the average over the spin_axis direction ; median_axis = the median over the spin_axis direction ; sigma = sigma for each spin period ; npoints = number of points in fit for each spin period ; sun_data = midpoint times of spitfit data ;keywords: ; plane_dim = Tells program which dimension to treat as the plane. 0=x, 1=y, 2=z. Default 0. ; axis_dim = Tells program which dimension contains axis to average over. Default 0. Will not ; create a tplot variable unless used with /spinaxis. ; min_points = Minimum number of points to fit. Default = 5. ; alpha = A parameter for finding fits. Points outside of sigma*(alpha + beta*i) ; will be thrown out. Default 1.4. ; beta = A parameter for finding fits. See above. Default = 0.4 ; phase_mask_starts = Time to start masking data. Default = 0 ; phase_mask_ends = Time to stop masking data. Default = -1 ; sun2sensor = Tells how much to rotate data to align with sun sensor. ; ;Example: ; thm_spinfit,'th?_fg?',/sigma ; 20-sep-2010, changed sign of sun2sensor to insure agressment between ; spinfit EFF and on-board EFS spinfit data, and between ; spinfit FGL and on-board FGS spinfit data. ;Written by Katherine Ramer ; $LastChangedBy: $ ; $LastChangedDate: $ ; $LastChangedRevision: $ ; $URL: $ ;- pro spinfit,arr_in_t,arr_in_data,arr_in_sunpulse_t,arr_in_sunpulse_data, $ A,B,C, avg_axis, median_axis,Sigma,Npoints,sun_data, $ min_points=min_points,alpha=alpha,beta=beta, $ plane_dim=plane_dim,axis_dim=axis_dim,phase_mask_starts=phase_mask_starts, $ phase_mask_ends=phase_mask_ends,sun2sensor=sun2sensor if not keyword_set(alpha) then alpha=1.4 else dprint, alpha if not keyword_set(beta) then beta = 0.4 if not keyword_set(min_points) then min_points=5 if not keyword_set(phase_mask_starts) then phase_mask_starts=0 if not keyword_set(phase_mask_ends) then phase_mask_ends=-1 if not keyword_set(sun2sensor) then sun2sensor=0 if not keyword_set(plane_dim) then plane_dim=0 if not keyword_set(axis_dim) then axis_dim=0 ; Make sure ARR_IN_T is monotonic (if necessary, remove 1 s chunks until it is). Update ARR_IN_DATA correspondingly: ; size_xxx=size(arr_in_t) monoton=1b non_monoton_detected = 0b k0=0L k1=k0+1L while monoton && ( k1 le size_xxx[1]-1 ) do begin if ( arr_in_t[ k1++ ] - arr_in_t[ k0++ ] lt 0 ) then begin non_monoton_detected = 1b monoton = 0b w = where( (arr_in_t ge (-0.5+arr_in_t[ k0-1 ]) ) and (arr_in_t le (0.5+arr_in_t[ k0-1 ]) ), complement = complement ) if complement[0] ne -1 then begin arr_in_t = temporary( arr_in_t[ complement ] ) arr_in_data = temporary( arr_in_data[ complement, *] ) endif size_xxx=size(arr_in_t) if w[0] ne -1 then begin if w[0] ne 0 then k0=w[0]-1 else k0 = 0 endif k1=k0+1L monoton=1b endif endwhile if non_monoton_detected then begin dprint, '*** WARNING: Non-monotonic time tags detected. 1 s chunks of data discarded until time tags are monotonic, but consult '+ $ 'THEMIS software team as to quality of time tags for these data!' endif ; find portion of data where the input overlaps the sunpulse times ; arr_in_sunpulse_t=arr_in_sunpulse_t+sun2sensor*arr_in_sunpulse_data/360 size_xxx=size(arr_in_t) overlap1=where(arr_in_sunpulse_t ge arr_in_t[0] and arr_in_sunpulse_t le arr_in_t[size_xxx[1]-1]) sizeoverlap=size(overlap1) dprint, "using axis dimension = ", axis_dim ; define dummy arrays to be filled later ; sigma=dblarr(sizeoverlap[1]-1) meany=dblarr(sizeoverlap[1]-1) result=dblarr(3) A=dblarr(sizeoverlap[1]-1) B=dblarr(sizeoverlap[1]-1) C=dblarr(sizeoverlap[1]-1) sun_data=dblarr(sizeoverlap[1]-1) avg_axis=dblarr(sizeoverlap[1]-1) median_axis=dblarr(sizeoverlap[1]-1) Npoints = ulonarr( sizeoverlap[1]-1 ) i=0L If(overlap1[0] Eq -1) Then Begin dprint, 'No good spin model' ;Bombs later Return Endif for i=0L, sizeoverlap[1]-2 do begin ; select a one period chunk of data: ; overlap = where(arr_in_t gt arr_in_sunpulse_t[overlap1[i]] and $ arr_in_t le arr_in_sunpulse_t[overlap1[i+1]], noverlap) if noverlap ge min_points then begin thx_xxx_keepy = arr_in_data[overlap, plane_dim] thx_xxx_keepz = arr_in_data[overlap, axis_dim] thx_xxx_keepx = arr_in_t[overlap] sun_data[i] = (arr_in_sunpulse_t[overlap1[i]]+arr_in_sunpulse_t[overlap1[i+1]])/2 ; make sure masked data is not selected ; sizemask = size(phase_mask_start) if max(phase_mask_ends)-min(phase_mask_starts) gt 360 then begin phase_mask_ends[0] = -1 dprint, "Warning: Phase Mask maxend-minstart > 360. No culling will occur." endif if not max(phase_mask_ends)-min(phase_mask_starts) lt 0 then begin phase = (2*!dPI/arr_in_sunpulse_data[overlap1[i]])*(thx_xxx_keepx[overlap] - arr_in_sunpulse_t[overlap1[i]]) for k = 0, sizemask[1]-1 do begin mask = where(phase gt phase_mask_starts[k] and phase lt phase_mask_ends[k]) thx_xxx_keepy[mask] = NaN thx_xxx_keepx[mask] = NaN thx_xxx_keepz[mask] = NaN endfor ;k endif ; throw out points 1.4 stddev (alpha) away from mean ; y = 0 meany[i] = mean(thx_xxx_keepy) sigma[i] = stddev(thx_xxx_keepy) nottoss = where(thx_xxx_keepy le meany[i]+sigma[i]*(alpha+beta*y) and thx_xxx_keepy ge meany[i]-sigma[i]*(alpha+beta*y)) If(nottoss[0] Ne -1) Then Begin sizenottoss = size(nottoss) sizelap = size(overlap) thx_xxx_keepy = thx_xxx_keepy[nottoss] thx_xxx_keepx = thx_xxx_keepx[nottoss] thx_xxx_keepz = thx_xxx_keepz[nottoss] ;print,sizelap[1]-sizenottoss[1], " initial points tossed" repeat begin ; recreate ranges after tossing out bad points above ; overlap = where(thx_xxx_keepx ge arr_in_sunpulse_t[overlap1[i]] and thx_xxx_keepx le arr_in_sunpulse_t[overlap1[i+1]]) ;Old location of performing least squares fit ; toss out bad points again and repeat if any points were tossed ; y = y+1 meany[i] = mean(thx_xxx_keepy) sigma[i] = stddev(thx_xxx_keepy) nottoss = where(thx_xxx_keepy le meany[i]+sigma[i]*(alpha+beta*y) and thx_xxx_keepy ge meany[i]-sigma[i]*(alpha+beta*y)) If(nottoss[0] Ne -1) Then Begin ; and thx_xxx_keepx ge arr_in_sunpulse_t[overlap1[i]] and thx_xxx_keepx le arr_in_sunpulse_t[overlap1[i+1]] ) thx_xxx_keepy = thx_xxx_keepy[nottoss] thx_xxx_keepx = thx_xxx_keepx[nottoss] thx_xxx_keepz = thx_xxx_keepz[nottoss] sizenottoss = size(nottoss) sizelap = size(overlap) NPoints[i] = sizenottoss[1] Endif Else Begin ;all points are bad for this spin? npoints[i] = 0 sizenottoss = size(nottoss) & sizenottoss[1] = 0 ;end this loop Endelse ;print,sizelap[1]-sizenottoss[1], " points tossed" endrep until (sizenottoss[1] eq sizelap[1] or sizenottoss[1] le min_points) avg_axis[i] = mean(thx_xxx_keepz) median_axis[i] = median(thx_xxx_keepz) Endif Else npoints[i] = 0 ;no good values if (Npoints[i] le min_points) then begin ; print, "No fit" A[i] = 'NaN' B[i] = 'NaN' C[i] = 'NaN' endif else begin ; calculate wt ; wt = (2*!dPI/arr_in_sunpulse_data[overlap1[i]])*(thx_xxx_keepx[overlap] - arr_in_sunpulse_t[overlap1[i]]) ; construct matrices ; sinwt = sin(wt) ;Repeated expressions. coswt = cos(wt) totsinwt = total(sinwt) totcoswt = total(coswt) sincoswt = sinwt*coswt totsincoswt = total(sincoswt) ; perform least squares fit ; result = la_least_squares([[sizelap[1], totcoswt, totsinwt], [totcoswt, total(coswt^2), totsincoswt], [totsinwt, totsincoswt, total(sinwt^2)]], $ transpose( [ total(thx_xxx_keepy), total(coswt*thx_xxx_keepy), total(sinwt*thx_xxx_keepy) ] ) ) ; print,result ; A[i] = result[0] B[i] = result[1] C[i] = result[2] endelse endif else begin A[i] = 'NaN' B[i] = 'NaN' C[i] = 'NaN' sun_data[i] = (arr_in_sunpulse_t[overlap1[i]]+arr_in_sunpulse_t[overlap1[i+1]])/2 endelse endfor ; i end