pro fit_plot_kcorr,$
    r=r,col_p=col_p,filter_title=filter_title,colour=colour,colmin=colmin,colmax=colmax,deg_z=deg_z,deg_col=deg_col,k_min=k_min,k_max=k_max,filename=filename,$
    good=good, kcorr_fit=bbb, kcorr_coeff=k1, plotrf=plotrf, col_rf=col_rf, c_rf_arr=c_rf_arr, c_rf_err=c_rf_err, c_rf_zmax=c_rf_zmax, c_rf_lab=c_rf_lab, c_rf_ls=c_rf_ls

if(n_elements(good) eq 0) then good=lindgen(n_elements(r))

tags = tag_names(r[0])
tag_col1 = where(tags eq 'AP3CORRMAG_'+strupcase(colour[0]), ctag_kcorr1)
tag_col2 = where(tags eq 'AP3CORRMAG_'+strupcase(colour[1]), ctag_kcorr2)
tag_kcorr = where(tags eq 'KCORR_'+strupcase(col_p), ctag_kcorr)
k1=fit_kcorr(r[good],col_p,colour=colour,deg_z=deg_z,deg_col=deg_col,kcfit=kcfit1,d_arr=d_arr1,good_sel=ggg1)
bbb=poly2d(r.z,r.(tag_col1)-r.(tag_col2),k1,deg1=deg_z,deg2=deg_col,/irreg)

if(n_elements(filename) ne 1) then filename='test'
if(strmid(filename,strlen(filename)-3) eq '.ps') then filename=strmid(filename,0,strlen(filename)-3)
if(strmid(filename,strlen(filename)-4) eq '.eps') then filename=strmid(filename,0,strlen(filename)-4)


m_k=median_line(r[good].z,r[good].(tag_kcorr)-bbb[good],min=0.03,max=0.6,bins=0.03,xnew=xx,yst=s_k)
set_plot,'ps'
!p.multi=[0,1,2]
device,filename=filename+'.ps',/color,/encaps,xs=9,ys=8
plot,[0,0.5999],[k_min,k_max],/nodata,xr=[0,0.5999],yr=[k_min,k_max],$
    xs=1,ys=1,ytitle='!5k!l'+col_p+'!n!3, mag',xtickn=[' ',' ',' ',' ',' ',' ',' '],ymar=[-4,1],xmar=[6,1]
;for i=0L,n_elements(bbb[good])-1 do oplot,[r[good[i]].z],[bbb[good[i]]],$
;      col=(bytscl(r[good[i]].(tag_col1)-r[good[i]].(tag_col2),colmin,colmax) < 254),psym=3
for i=0L,n_elements(bbb[good])-1 do oplot,[r[good[i]].z],[r[good[i]].(tag_kcorr)],$
      col=(bytscl(r[good[i]].(tag_col1)-r[good[i]].(tag_col2),colmin,colmax) < 254),psym=3

if(keyword_set(plotrf) and (n_elements(col_rf) eq n_elements(r))) then begin
    if(n_elements(c_rf_arr) lt 1) then c_rf_arr=[0.77]
    if(n_elements(c_rf_err) ne n_elements(c_rf_arr)) then c_rf_err=c_rf_arr*0.0+0.04
    if(n_elements(c_rf_zmax) ne n_elements(c_rf_arr)) then c_rf_zmax=c_rf_arr*0.0+0.6
    if(n_elements(c_rf_ls) ne n_elements(c_rf_arr)) then c_rf_ls=c_rf_arr*0+2

    for k=0,n_elements(c_rf_err)-1 do begin
        gal_sel = where(abs(col_rf - c_rf_arr[k]) le c_rf_err[k], cgal_sel)
        if(cgal_sel gt 100) then begin
            m_l_gal=median_line(r[gal_sel].z,r[gal_sel].(tag_kcorr),min=0.03,max=c_rf_zmax[k],bins=0.03,xnew=xx1,yst=s_m_l_gal)
            oplot,xx1,m_l_gal,linest=c_rf_ls[k]
            if(n_elements(c_rf_lab) gt k) then begin
                xyouts,xx1[n_elements(xx1)-1]+0.0,m_l_gal[n_elements(xx1)-1],c_rf_lab[k],align=0.0
            endif
        endif
    endfor
endif

plot,xr=[0,0.5999],yr=[-0.1499,0.1499],xx,m_k,xtitle='Redshift',ytitle='!5k!l!3comp!n - !5k!l!3fit!n',ymar=[3,4],xs=1,ys=1,xmar=[6,1],$
    ytickn=['',' ','',' ','']
oplot,[0,0.6],[0,0],linest=0,col=40
oplot,xx,m_k-s_k,linest=1
oplot,xx,m_k+s_k,linest=1
colbar=findgen(254)
tv,colbar,2.0,6.7,xsize=1.3,ys=0.2,/cent
plot,[colmin,colmax],[0,1],/nodata,/noer,ticklen=0,position=1000.0*[2.0,6.7,2.0+1.3,6.7+0.2],$
    /dev,xs=1,ys=1,ytickn=[' ',' ',' ',' ',' ',' ',' ',' '],xticks=1,chars=0.7,yticks=1
xyouts,1000.0*2.65,1000.0*7.0,'!5'+colour[0]+' - '+colour[1]+'!3',/dev,align=0.5
xyouts,1000.0*1.5,1000.0*7.3,filter_title,/dev
device,/close
!p.multi=0
set_plot,'x'

deg1=deg_z
deg2=deg_col
c=k1
s=size(c)
if(s[0] eq 2) then begin
    deg1=s[1]
    deg2=s[2]
    c1=c
endif else begin
    if((n_elements(deg1) ne 1) and (n_elements(deg2) ne 1)) then deg2=3 ;;;; default, horribly wrong :)
    deg1 = n_elements(c)/(float(deg2) + 1) + float(deg2)/2 - 1
    c1=dblarr(deg1+1,deg2+1)
    ii=0L
    for i=0, deg1 do begin
        for j=0, deg2 do begin
            if(i+j gt (deg1>deg2)) then continue
            c1[i,j]=c[ii]
            ii=ii+1L
        endfor
    endfor
endelse

writefits,filename+'.fits',c
h=headfits(filename+'.fits')
sxaddpar,h,'SERVICE','K correction 2D polynomial approximation'
sxaddpar,h,'DEG_Z',deg_z
sxaddpar,h,'DEG_COL',deg_col
sxaddpar,h,'BAND',col_p
sxaddpar,h,'COLOUR',colour[0]+'-'+colour[1]
writefits,filename+'.fits',c,h

col_title=colour[0]+' - '+colour[1]
openw,u,filename+'.txt',/get_lun
printf,u,'\begin{table}'
;printf,u,'\caption{Coefficients of the two-dimensional polynomial approximation of $k$-corrections for the $'+filter_title+'$ band.\label{tab'+colour[0]+'_'+colour[1]+'}}'
printf,u,'\caption{Same as in Table~A1 but for the '+filter_title+' band as a function of redshift and the $'+col_title+'$ colour.\label{tab'+colour[0]+'_'+colour[1]+'}}'
printf,u,'\begin{tabular}{l|cccc}'
printf,u,' & $('+col_title+')^0$ & $('+col_title+')^1$ & $('+col_title+')^2$ & $('+col_title+')^3$ \\'
printf,u,'\hline'
for i=0,deg_z do begin
  printf,u,'$Z^{'+string(i,format='(i1)')+'}$ & ',c1[i,0],' & ',c1[i,1],' & ',c1[i,2],' & ',c1[i,3],'\\',$
        format='(%"%s %g %s %g %s %g %s %g %s")'
endfor
printf,u,'\hline'
printf,u,'\end{tabular}'
printf,u,'\end{table}'
close,u
free_lun,u


end