;20180404 Ali
;calculates the fraction of the sep fov covered by mars surface and mars shine
;rmars: mars radius (km)
;posmar: center of mars position w.r.t maven (km)
;possun: unit vector pointing to the sun from maven (or mars)
;fov: if set, calculations are done for only within the fov of each sep
;resdeg: angular resolution of calculations in degrees
;response: if set, plots fov response
;vector: if set, vectorization is used (much faster, but much higher memory usage)
;
function mvn_sep_fov_mars_shine,rmars,posmar,possun,fov=fov,resdeg=resdeg,response=response,vector=vector

  fnan=!values.f_nan
  srefa=15.5 ;ref (phi) angle
  scrsa=21.0 ;cross (theta) angle

  if ~keyword_set(resdeg) then resdeg=1.
  phid=resdeg*findgen(360./resdeg)-180. ;degrees
  thed=resdeg*findgen(180./resdeg)
  spc=[-45.,+45.,+135.,-135.] ;SEP 1F,2F,1R,2R phi center
  wthe=where(abs(thed-90.) le scrsa,/null,nth2)
  if keyword_set(fov) then begin
    wphi=where((abs(phid-spc[0]) le srefa) or (abs(phid-spc[1]) le srefa) or (abs(phid-spc[2]) le srefa) or (abs(phid-spc[3]) le srefa),/null)
    phid=phid[wphi]
    thed=thed[wthe]
    wthe=where(abs(thed-90.) le scrsa,/null,nth2)
  endif
  nphi=n_elements(phid)
  nthe=n_elements(thed)
  fthe=(cos(!dtor*(thed[wthe]-90.))-cos(!dtor*scrsa))/(1.-cos(!dtor*scrsa)) ;weighting based on detector angular response

  x=cos(!dtor*(phid+45.))#sin(!dtor*thed)
  y=sin(!dtor*(phid+45.))#sin(!dtor*thed) ;sep1z
  z=replicate(1.,nphi)#cos(!dtor*thed) ;-sep1y

  sizemar=size(posmar,/dim)
  if n_elements(sizemar) eq 1 then nt=1 else nt=sizemar[1]
  mars_surfa=replicate(fnan,[4,nt])
  mars_shine=mars_surfa
  atmo_shine=mars_surfa
  mshine_fov=mars_surfa
  ashine_fov=mars_surfa

  if keyword_set(vector) then begin ;vectorization: faster, but uses more memory!
    dvec=[transpose(reform(x,nphi*nthe)),transpose(reform(z,nphi*nthe)),transpose(reform(y,nphi*nthe))]
    mvn_sep_fov_mars_incidence,rmars,posmar,dvec,spoint
    cossza=reform(total((spoint-transpose(rebin(posmar,[3,nt,nphi*nthe]),[0,2,1]))*transpose(rebin(possun,[3,nt,nphi*nthe]),[0,2,1]),1),[nphi,nthe,nt])/rmars
    tanalt=cossza
  endif else begin
    cossza=replicate(fnan,[nphi,nthe,nt])
    tanalt=cossza
    atmosh=cossza
    radmar=sqrt(total(posmar^2,1)) ;radial distance of MAVEN from Mars (km)
    for ithe=0,nthe-1 do begin
      for iphi=0,nphi-1 do begin
        dvec=[x[iphi,ithe],z[iphi,ithe],y[iphi,ithe]]
        pdmmar=total(posmar*rebin(dvec,[3,nt]),1)/radmar
        talt=radmar*sqrt(1.d0-pdmmar^2)-rmars ;tangent altitude of dvec (km)
        wtalt0=where(talt lt 0. and pdmmar gt 0.,/null,nws) ;surface
        wtagt0=where(talt gt 0. and talt lt 100. and pdmmar gt 0.,/null,nwt) ;for atmosphere shine
        talt[where(pdmmar lt 0. or talt lt 0.,/null)]=fnan  ;where dvec away from Mars or negative tanalt
        tanalt[iphi,ithe,*]=talt
        if nws gt 0 then begin
          mvn_sep_fov_mars_incidence,rmars,posmar[*,wtalt0],dvec,spoint
          cossza[iphi,ithe,wtalt0]=total((reform(spoint)-posmar[*,wtalt0])*possun[*,wtalt0],1)/rmars
        endif
;        if nwt gt 0 then begin
          postal=dvec#(radmar*pdmmar)-posmar ;tangent altitude vector from Mars center (km)
          cosszatal=total(postal*possun,1)/sqrt(total(postal^2,1)) ;cos(sza of tanalt)
          atmosh[iphi,ithe,*]=(cosszatal gt 0.)*(10.+cosszatal)*exp(-talt/10.)
;        endif
      endfor
    endfor
  endelse
  for isep=0,3 do begin
    wph2=where(abs(phid-spc[isep]) le srefa,/null,nph2)
    fphi=(cos(!dtor*(phid[wph2]-spc[isep]))-cos(!dtor*srefa))/(1.-cos(!dtor*srefa)) ;weighting based on detector angular response
    wfov=replicate(1.,[nph2,nth2]) ;constant (uniform) weighting for fov elements
    wfo2=fphi#fthe
    if keyword_set(response) then p=image(wfo2,phid[wph2],thed[wthe]-90.,rgb=33,/o,min=0.,max=1.)
    cosszafov=cossza[wph2,wthe,*] ;elements within fov
    atmoshfov=atmosh[wph2,wthe,*]
    mars_surfa[isep,*]=total(total(finite(cosszafov)*rebin(wfov,[nph2,nth2,nt]),1),1)/total(wfov) ;mars surface (disc)
    mars_shine[isep,*]=total(total((cosszafov gt 0.)*cosszafov*rebin(wfov,[nph2,nth2,nt]),/nan,1),1)/total(wfov) ;mars shine
    atmo_shine[isep,*]=total(total((atmoshfov gt 0.)*atmoshfov*rebin(wfov,[nph2,nth2,nt]),/nan,1),1)/total(wfov) ;atmo shine
    mshine_fov[isep,*]=total(total((cosszafov gt 0.)*cosszafov*rebin(wfo2,[nph2,nth2,nt]),/nan,1),1)/total(wfo2) ;mars shine convolved w/ fov response
    ashine_fov[isep,*]=total(total((atmoshfov gt 0.)*atmoshfov*rebin(wfo2,[nph2,nth2,nt]),/nan,1),1)/total(wfo2) ;mars shine convolved w/ fov response
  endfor

  fraction={fov:['SEP1F','SEP2F','SEP1R','SEP2R'],mars_surfa:mars_surfa,mars_shine:mars_shine,atmo_shine:atmo_shine,mshine_fov:mshine_fov,ashine_fov:ashine_fov,cossza:cossza,tanalt:tanalt,atmosh:atmosh,phid:phid,thed:thed}
  return,fraction

  if 0 then begin ;old method: Mars shine on surface based on cos(sza)
    x=sin(th/2.)#cos(th)
    y=sin(th/2.)#sin(th)
    z=rebin(cos(th/2.),[nth,nth])
    xyz=reform(transpose([[[x]],[[y]],[[z]]]),[3,nth^2])
    cossza=total(xyz*rebin(possun,[3,nth^2]),1)
    xyzmar=rmars*xyz+rebin(posmar,[3,nth^2])
    xyzdir=xyzmar/(ones3#sqrt(total(xyzmar^2,1)))
    dist=sqrt(total(xyzmar^2,1)) ;distance of point on the mars surface from maven
    wdlt=where(dist lt sqrt(mvnrad^2-rmars^2),/null)

    septheta=90.-!radeg*acos(xyzdir[2,wdlt]) ;sep-xy angle (degrees)
    sep1fphi=!radeg*atan(xyzdir[1,wdlt],xyzdir[0,wdlt]) ;sep1f-xz angle (degrees) [-180,180]
    sepphi=sep1fphi-45. ;to align phi=0 with s/c +Z axis
    wlt180=where(sepphi lt -180.,/null)
    if n_elements(wlt180) gt 0 then sepphi[wlt180]+=360.
    fovmars=replicate(0.,[360,180])
    fovshin=fovmars
    ph=floor(sepphi)+180  ;[0,359] ;ref angle
    th=floor(septheta)+90 ;[0,179] ;cross angle
    shine=cossza
    shine[where(shine lt 0.,/null)]=0.
    np=n_elements(cossza) ;number of points
    for ip=0,np-1 do begin
      fovmars[ph[ip],th[ip]]+=1.
      fovshin[ph[ip],th[ip]]+=shine[ip]
    endfor
    marshin=fovshin/fovmars
    p=image(marshin,findgen(360)-180.,findgen(180)-90.,rgb=colortable(64,/reverse),/o,min=0.,max=1.)
    szascaled=bytscl(1.-cossza,min=0.,max=1.)
    p=scatterplot(/o,sepphi,septheta,rgb=64,sym='o',/sym_filled,magnitude=szascaled)
  endif
end