;
; Multi guidestar Tip/Tilt error calculator
;
;  Given a number of tip/tilt stars calculate:
;   1. From the brightness, wavelength, and exposure time compute the SNR
;      and appropriate noise-equivalent-angles
;   2. Using the star positions relative to a science object compute the
;      cross-correlation factors, optimal gain matrix, and rms residual
;      tip/tilt at the science object
;
pro TipTiltError,Cn2,h,D,ttStarList,scienceStarList,coef,tiltNoise,G,P
	nstars = (size(ttStarList))[2]
;
;  calculate the tilt variance given D and Cn2
	tilt_variance = 3.03801*Cn2/d^(1./3.) ; radians squared "signal" from atmosphere
;
;  cross-correlations amongst tip/tilt stars
;
	Rxx = fltarr(nstars*2,nstars*2)
	Rxy = fltarr(2,nstars*2)
	cxx = Rxx
	cxy = Rxy

	nlayers = (size(h))[1]
	for layer = 0,nlayers-1 do begin
		for star1 = 0,nstars-1 do begin
		  index1 = 2*star1
		  for star2 = 0,nstars-1 do begin
		    index2 = 2*star2
		    s = (ttStarList[*,star1]-ttStarList[*,star2])*h[layer]/(D/2)
			sn = norm(s)
			if (sn ne 0) then begin
			  cth = s[0]/sn
			  sth = s[1]/sn
			endif else begin
			  cth = 1.0
			  sth = 0.0
			endelse
			cx = interpol(coef[1,*],coef[0,*],[sn])
			cy = interpol(coef[2,*],coef[0,*],[sn])
			cc = diagonal([cx,cy])
			rm = [[cth,sth],[-sth,cth]]
		    cc = rm ## cc ## transpose(rm)
		    cxx[index1:index1+1,index2:index2+1] = cc
		  endfor
	endfor
;
; cross-correlations of tip/tilt stars to science object at [scienceStarLoc]
;
	for star = 0,nstars-1 do begin
	  index = 2*star
	  s = (ttStarList[*,star]-scienceStarList)*h[layer]/(D/2)
	  sn = norm(s)
	  if (sn ne 0) then begin
	    cth = s[0]/sn
	    sth = s[1]/sn
	  endif else begin
	    cth = 1.0
	    sth = 0.0
	  endelse
	  cx = interpol(coef[1,*],coef[0,*],[sn])
	  cy = interpol(coef[2,*],coef[0,*],[sn])
	  cc = diagonal([cx,cy])
	  rm = [[cth,sth],[-sth,cth]]
	  cc = rm ## cc ## transpose(rm)
	  cxy[*,index:index+1] = cc
	endfor

	;
	; create covariance matrices Rxx=variance*cxx+noiseVariance and Ryx=variance*cyx
	;
	Rxx += cxx*tilt_variance[layer]
	Rxy += cxy*tilt_variance[layer]

	endfor ; layers

	for star = 0,nstars-1 do begin
	  index = 2*star
	  Rxx[index:index+1,index:index+1] += diagonal([1,1]*tiltNoise[star]^2)
	endfor
	;
	; compute gain matrix and error covariance matrix
	;
	total_tilt_variance = total(tilt_variance)
	Rxx /= total_tilt_variance
	Rxy /= total_tilt_variance
	G = transpose(Rxy) ## invert(Rxx)
	P = Identity(2) - G ## Rxy
;
        return
end
; -------------------------------------------------------------------------------------
microns = 1.e-6
arcsec = !pi/(180.*3600.)
arcmin = 60.*arcsec
degrees = 60.*arcmin
km = 1.e3
; ---------------------------
; cases:
;field_set = ['goods_north','goods_south','cosmos','goods_north','goods_south','cosmos','goods_north','goods_south','cosmos']
;outerScale_set = [0,0,0,2,2,2,7,7,7]
;caseText_set = ['goodsN_0','goodsS_0','cosmos_0','goodsN_20','goodsS_20','cosmos_20','goodsN_70','goodsS_70','cosmos_70']
field_set = ['goods_south']
outerScale_set = [0]
caseText_set = ['goodsS_0']
ncases = (size(field_set))[1]
for theCase = 0,ncases-1 do begin
	caseText = caseText_set[theCase]
	outerScale = outerScale_set[theCase]
	field = field_set[theCase]
;caseText = 'goodsS_20'
;  -- problem parameters --
D = 10.
r0 = 0.2
lambda0 = 0.55*microns
h = [0,1.8,3.3,5.8,7.4,13.1,15.8]*km
Cn2_profile = [0.646,0.078,0.119,0.035,0.025,0.082,0.015]; [0,0,0,0,0,0,1];
C0 = r0^(-5./3.)*(lambda0/(2*!pi))^2/0.423
Cn2 = C0*Cn2_profile
;field = 'goods_south';'goods_south';'goods_north';'cosmos'
fakeStars = 0
nstars = 24 ; (fake stars only)
fakeStarFieldRadius = 2*arcmin
n_ttStars = 5  ; number of tip/tilt stars to use
;outerScale = 2 ; number of telescope diameters, [0=infinity, 1, 2, 7]
ttBand = 'J'
ttStrehl = 0.5
exposureTime = 0.01
writeFiles = 0
; ---------------------------
calcBandData,bandName,bandData
ind = (where(ttBand eq bandName))[0]
lambdaTT = bandData[0,ind]*microns
airyCore = lambdaTT/D
seeingFWHM = lambdaTT/(r0*(lambdaTT/lambda0)^(6./5.))
;clearwins
;
;  -- read in tilt correlation coefficients --
;
dir = 'C:\Documents and Settings\Don Gavel\My Documents\Adaptive Optics\Keck\Keck NGAO\Keck NGAO Extragalactic Science Modeling\'
file = 'tipTiltCorrelationCoefficients.dat'
u = read_ascii(dir+file,comment_symbol=';')
u = u.field01
plotCC = 0
if (plotCC) then begin
	window,0
	plot,u[0,*],u[3,*],xrange=[0,4],yrange=[-1,1]
	oplot,u[0,*],u[4,*],linestyle=2
	oplot,u[0,*],u[5,*]
	oplot,u[0,*],u[6,*],linestyle=2
	oplot,u[0,*],u[7,*]
	oplot,u[0,*],u[8,*],linestyle=2
	oplot,u[0,*],u[9,*]
	oplot,u[0,*],u[10,*],linestyle=2
endif
k = 3 ; by default, infinite outer scale
alpha = 3.039
if (outerScale eq 1) then begin
	k = 9
	alpha = 0.0537852
endif
if (outerScale eq 2) then begin
	k = 7
	alpha = 0.209782
endif
if (outerScale eq 7) then begin
	k = 5
	alpha = 0.823902
endif
ttCoef = [u[0,*],u[k,*],u[k+1,*]]
;
;  -- calculate the initial rms tip/tilt --
;
tilt_variance = alpha*Cn2/D^(1./3.) ; radians squared "signal" from atmosphere
total_tilt_variance = total(tilt_variance)
;
if (fakeStars) then begin
	;  -- generate a fake star list --
	starLocs = (randomu(seed,2,nstars)-0.5)*2*fakeStarFieldRadius
	tiltNoise = 0.001*(ones(nstars))*arcsec ; radians rms measurement noise
endif
if (not fakeStars) then begin
	;
	;  -- read in the field stars
	;
	dir = 'C:\Documents and Settings\Don Gavel\My Documents\Adaptive Optics\Keck\Keck NGAO\Keck NGAO Extragalactic Science Modeling\'
	if (field eq 'goods_north') then begin
		file = 'GoodsStarLists_fromJasonMelborne\goods_n_lgs.dat'
		s = read_ascii(dir+file,comment_symbol=';')
		s = s.field1
		nstars = (size(s))[2]
		starField = replicate('GOOD',nstars)
		ra = reform(s[1,*]+s[2,*]/60.+s[3,*]/3600.)
		dec = reform(sign(s[4,*])*(abs(s[4,*])+s[5,*]/60.+s[6,*]/3600.))
		ra0 = max(ra)
		dec0 = min(dec)
		bmag = s[7,*]
		rmag = s[8,*]
	endif
	if (field eq 'goods_south') then begin
		file = 'GoodsStarLists_fromJasonMelborne\goods_s_gems.dat'
		template = {version:1.0, datastart:long(0), delimiter:byte(32), missingvalue:0.0, commentsymbol:';',fieldcount:long(10), $
					fieldtypes:[4,4,4,4,4,4,4,4,4,7],fieldnames:['FIELD01','FIELD02','FIELD03','FIELD04','FIELD05','FIELD06','FIELD07','FIELD08','FIELD09','FIELD10'], $
					fieldlocations:[4, 10, 13, 16, 27, 31, 34, 41, 47, 52], fieldgroups:[0,0,0,0,0,0,0,0,0,9]}
		s = read_ascii(dir+file,template=template)
		starField = s.field10
		s = s.field01
		ra = reform(s[1,*]+s[2,*]/60.+s[3,*]/3600.)
		dec = reform(sign(s[4,*])*(abs(s[4,*])+s[5,*]/60.+s[6,*]/3600.))
		ra0 = max(ra)
		dec0 = min(dec)
		bmag = s[7,*]
		vmag = s[8,*]
	endif
	if (field eq 'cosmos') then begin
		file = 'usno-search-COSMOS.txt'
		u = read_ascii(dir+file,comment=';')
		u = u.field01
		ra = reform(u[8,*])/15. ; this is actually RA converted to degrees
		dec = reform(u[9,*]) ; degrees
		nstars = (size(ra))[1]
		mag1 = fltarr(nstars)
		band1 = strarr(nstars)
		mag2 = fltarr(nstars)
		band2 = strarr(nstars)
		bands = ['B','R','B','R','I']
		b1mag = u[16,*] ; B band magnitude
		r1mag = u[17,*] ; R band magnitude
		b2mag = u[18,*]
		r2mag = u[19,*]
		ind = where(b1mag eq 99.99,count)
		if (count ne 0) then b1mag[ind] = b2mag[ind]
		ind = where(r1mag eq 99.99,count)
		if (count ne 0) then r1mag[ind] = r2mag[ind]
		rmag = r1mag
		bmag = b1mag
		imag = u[20,*] ; I band magnitude
		goodStar = ones(nstars)
		ind = where(rmag eq 99.99,count)
		if (count ne 0) then goodStar[ind] = 0
		ind = where(bmag eq 99.99,count)
		if (count ne 0) then goodStar[ind] = 0
		ind = where(imag eq 99.99,count)
		if (count ne 0) then goodStar[ind] = 0
		; -- trim down the number of stars to a reasonable number in a sub field --
		;  (this will end up being a statistical sample)
		cx = 150.085/15.
		cy = 2.27
		fovx = (150.147-150.025)/15.
		fovy = 2.33-2.22
		u = where( goodStar and (ra le (cx+fovx/2.)) and (ra ge (cx-fovx/2.)) and (dec le (cy+fovy/2.)) and (dec ge (cy-fovy/2.)) )
		ra = ra[u]
		dec = dec[u]
		bmag = bmag[u]
		rmag = rmag[u]
		imag = imag[u]
		ra0 = max(ra)
		dec0 = min(dec)
		nstars = (size(ra))[1]
		starField = replicate('COSMOS',nstars)
	endif
	nstars = (size(ra))[1]
	starLocs = fltarr(2,nstars)
	starFlux = fltarr(nstars)
	starTemp = fltarr(nstars)
	tiltNoise = fltarr(nstars)
	cdec = -cos(dec0*degrees)
	for star = 0,nstars-1 do begin
		starLocs[0,star] = ra[star]*15*cdec; *degrees
		starLocs[1,star] = dec[star]; *degrees
		if (field eq 'goods_north') then begin
			starTemp[star] = photoColor('B',bmag[star],'R',rmag[star]) ; K
			starFlux[star] = fluxRatio(starTemp[star],'J','R',/photons)*photoCounts('R',rmag[star]) ; ph/s/m2
		endif
		if (field eq 'goods_south') then begin
			starTemp[star] = photoColor('B',bmag[star],'V',vmag[star]) ; K
			starFlux[star] = fluxRatio(starTemp[star],'J','V',/photons)*photoCounts('V',vmag[star]) ; ph/s/m2
		endif
		if (field eq 'cosmos') then begin
			starTemp[star] = photoColor('B',bmag[star],'R',rmag[star]) ; K
			starFlux[star] = fluxratio(starTemp[star],'J','R',/photons)*photoCounts('R',rmag[star]) ; ph/s/m2
		endif
		tiltNoise[star] = airyCore/sqrt(ttStrehl*starFlux[star]*exposureTime*D^2*(!pi/4.))
	endfor
;	tiltNoise = 0.001*(ones(nstars))*arcsec ; radians rms measurement noise
endif

; -- plot the stars in the star list --

cx = average(starLocs[0,*])
cy = average(starLocs[1,*])
fovx = max(starLocs[0,*])-min(starLocs[0,*])
fovy = max(starLocs[1,*])-min(starLocs[1,*])
fov = max([fovx,fovy])

;fov = max(abs(starLocs))*1.1
doPlot = 1
if (doPlot) then begin
	if (field eq 'goods_north') then win = 0
	if (field eq 'goods_south') then win = 1
	if (field eq 'cosmos') then win = 2
    window,win,xsize=500,ysize=500,title = field
;    plot,starLocs[0,*]/arcmin,starLocs[1,*]/arcmin,xtitle='arcmin',ytitle='arcmin',psym=2,xrange=([-1,1]*fov/2+cx)/arcmin,yrange=([-1,1]*fov/2+cy)/arcmin
	plot,starLocs[0,*]/cdec,starlocs[1,*],/nodata,xtitle='RA',ytitle='DEC',xrange=([-1,1]*fov/2+cx)/cdec,yrange=([-1,1]*fov/2+cy)
	for star = 0,nstars-1 do begin
		if (starField[star] eq 'GOOD') then oplot,[starLocs[0,star]]/cdec,[starLocs[1,star]],psym=1
		if (starField[star] eq 'GEMS') then oplot,[starLocs[0,star]]/cdec,[starLocs[1,star]],psym=4
		if (starField[star] eq 'COSMOS') then oplot,[starLocs[0,star]]/cdec,[starLocs[1,star]],psym=1
	endfor

;    for star = 0,nstars-1 do begin
;        xyouts,starLocs[0,star]/arcmin+.05,starLocs[1,star]/arcmin,strtrim(string(star),2)
;    endfor
endif
;
; --- sample the field uniformly and compute the tip/tilt residual at each field point ---
;
nFieldAcross = 21 ; pick an odd number
n2 = (nFieldAcross-1)/2
nFieldPos = nFieldAcross^2
fieldLocsX = (ones(nFieldAcross) ## findgen(nFieldAcross) - n2)/n2*fov*0.8
fieldLocsY = transpose(fieldLocsX)
fieldLocsX += cx
fieldLocsY += cy
scienceStarLoc = fltarr(2,1)
corrected_tilt_sigma = fltarr(nFieldAcross,nFieldAcross)
print,'uncorrected  tilt sigma ',sqrt(total_tilt_variance)/arcsec,' arcsec rms'; rss sum tip/tilt over Cn2
print,'tilt noises ',tiltNoise/arcsec/.001,' milli-arcsec rms'
print,'shears ',h*fov/D,' metapupil diameters'
for fieldPos = 0,nFieldPos-1 do begin
;
; -- rank order the tip/tilt stars by distance and pick the closest n_ttStars ---
;
    scienceStarLoc = [fieldLocsX[fieldPos],fieldLocsY[fieldPos]]
    starDist = starLocs - (ones(nstars) ## scienceStarLoc)
    starDist = sqrt((starDist*starDist) ## ones(2))
    index = sort(starDist)
    nearStarLocs = starLocs[*,index[0:n_ttStars-1]]
    if (doPlot) then oplot,[scienceStarLoc[0]]/cdec,[scienceStarLoc[1]],psym=3,color=100
    TipTiltError,Cn2,h,D,nearStarLocs*degrees,scienceStarLoc*degrees,ttCoef,tiltNoise,G,P
;
; -- report results --
;
    corrected_tilt_sigma[fieldPos] = norm(sqrt(diagonal(p)))*sqrt(total_tilt_variance)
    print,'corrected tilt sigma ',norm(sqrt(diagonal(p))*sqrt(total_tilt_variance)/arcsec/.001),' milli-arsec rms'
endfor
fx0 = min(fieldLocsX)/cdec
fy0 = min(fieldLocsY)
dfx = (fieldLocsX[1,0]-fieldLocsX[0,0])/cdec
dfy = fieldLocsY[0,1]-fieldLocsY[0,0]
disp,corrected_tilt_sigma/arcsec/.001,caseText,x0=[fx0,fy0],dx=[dfx,dfy]

outputFile = caseText+'.fits'
outputFolder = 'ttErrorMaps\'
mkhdr,outHeader,corrected_tilt_sigma
sxaddpar,outHeader,'UNITS','milli-arcseconds tip/tilt error'
sxaddpar,outHeader,'INSTR','Keck NGAO'
sxaddpar,outHeader,'FIELD',field,'star field'
sxaddpar,outHeader,'RA_0',fx0,'lower left pixel RA, degrees'
sxaddpar,outHeader,'DEC_0',fy0,'lower left pixel DEC, degrees'
sxaddpar,outHeader,'DRA',dfx,'delta RA, degrees'
sxaddpar,outHeader,'DDEC',dfy,'delta DEC, degrees'
sxaddpar,outHeader,'FILTER','J','TT star sensing band'
sxaddpar,outHeader,'EXPTIME',exposureTime,'exposure time, seconds'
sxaddpar,outHeader,'STREHL',ttStrehl,'TT star Strehl'
sxaddpar,outHeader,'OUTSCL',outerScale*D,'outer scale, meters'

if (writeFiles) then writefits,dir+outputFolder+outputFile,corrected_tilt_sigma/arcsec/.001,outHeader


endfor ; for ... ncases

end