pro spinfit,arr_in_t,arr_in_data,arr_in_sunpulse_t,arr_in_sunpulse_data, $
           A,B,C,spin_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 print,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



;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)

print,"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)

for i=0L,n_elements(overlap1)-2 do begin
;select a one period chunk of data
overlap=where(arr_in_t ge arr_in_sunpulse_t[overlap1[i]] and arr_in_t le arr_in_sunpulse_t[overlap1[i+1]])

if n_elements(overlap) 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
  print,"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
endif

;throw out points 1.4 stddev 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))
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]])
x=size(overlap)

;calculate wt
wt=(2*!dPI/arr_in_sunpulse_data[overlap1[i]])*(thx_xxx_keepx[overlap] - arr_in_sunpulse_t[overlap1[i]])

;construct matrices
row1=[x[1],total(cos(wt)),total(sin(wt))]
row2=[total(cos(wt)),total(cos(wt)*cos(wt)),total(sin(wt)*cos(wt))]
row3=[total(sin(wt)),total(sin(wt)*cos(wt)),total(sin(wt)*sin(wt))]
matrix=[[row1],[row2],[row3]]
bmatrix=[total(thx_xxx_keepy),total(cos(wt)*thx_xxx_keepy),total(sin(wt)*thx_xxx_keepy)]
bmatrix=transpose(bmatrix)

;perform least squares fit

result=la_least_squares(matrix,bmatrix)
;print,result
A[i]=result[0]
B[i]=result[1]
C[i]=result[2]


;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))
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=sizenottoss[1]
;print,sizelap[1]-sizenottoss[1], " points tossed"

endrep until (sizenottoss[1] eq sizelap[1] or Npoints le min_points)

avg_axis[i]=mean(thx_xxx_keepz)
median_axis[i]=median(thx_xxx_keepz)

if (Npoints le min_points) then begin
print,"No fit"
A[i]='NaN'
B[i]='NaN'
C[i]='NaN'
endif
endif else begin
A[i]='NaN'
B[i]='NaN'
C[i]='NaN'
endelse
endfor
end