;
; PURPOSE: Produce data file for harmonic coefficients using photospheric
;          magnetic fields
; NOTE: The obtained coefficients can be used to both PFSS and HCCSSS models
; WRITTEN: Xuepu Zhao Jan. 2006
;
pro ghf_phof06,phof,nmax,ga,ha,ghf,WSO=wso

    if N_params() LT 1 then begin
       print,'pro ghf_phof06,phof,nmax,ga,ha,ghf'
       return
    endif
    if keyword_set(wso) then bra=readwsosc(phof,/cali,/br) $
                        else bra=readfits(phof,hd)
    sz=SIZE(bra)
    if keyword_set(wso) then xsz=sz(1)-1 else xsz=sz(1) 
    ysz=sz(2)
    zgrid06,xsz,ysz,phd,thd,lad,cth,sth
    pdp06,cth,nmax,pp,dp
    csmph06,phd,nmax,cmph,smph
;
    get_lun,unit
    openw,unit,ghf
    printf,unit,phof
    printf,unit,nmax
;
    ga=fltarr(nmax+1,nmax+1)
    ha=ga
    ij=xsz*ysz
    for n=0,nmax do begin         ; n=0 monopole at photosphere 25Jan2006
        coe=(2*n+1.0)/ij
        for m=0,n do begin
            ga(n,m)=0
            ha(n,m)=0
            for i=0,xsz-1 do begin
            ;   if i gt 50 then print,i
                for j=0,ysz-1 do begin
                    ga(n,m)=ga(n,m)+bra(i,j)*pp(n,m,j)*cmph(m,i)
                    ha(n,m)=ha(n,m)+bra(i,j)*pp(n,m,j)*smph(m,i)
                endfor
            endfor
            ga(n,m)=coe*ga(n,m)
            ha(n,m)=coe*ha(n,m)
            printf,unit,n,m,ga(n,m),ha(n,m)
        endfor
    endfor
    free_lun,unit

end