;
pro test_w05sc,ps=ps
common saveSCplm,prevth0,ColatTable,PlmTable,nlms
;
 prevth0 = 0. ; init
;
; Compile these 2 source files from command line before
; calling this test driver:
;
; .r w05sc.pro
; .r w05scDerived.pro
;
 make_plots = 1
;
 mlat = [-90., -88.1238292398491, -86.2386359278657, -84.3344382773342,$
    -82.4013318763435, -80.4295344892688, -78.4094552099168,$
    -76.331796630125, -74.1876988925388, -71.9689341802758,$
    -69.6681589022773, -67.2792279882741, -64.7975706790533,$
    -62.2206194320588, -59.5482728298363, -56.7833601290164,$
    -53.9320608459732, -51.0042204168578, -48.0134966005524,$
    -44.9772754602266, -41.916313892128, -38.8540980954293,$
    -35.8159497801506, -32.8279553674349, -29.9158266703621,$
    -27.1038148776609, -24.4137889090065, -21.8645574169981,$
    -19.4714697638694, -17.2462861630082, -15.1972697734841,$
    -13.3294282264571, -11.6448185129562, -10.142824406667,$
    -8.82031765103987, -7.67162666281269, -6.68827297583048,$
    -5.85851734698832, -5.16689314460211, -4.5940469432968,$
    -4.11722526306697, -3.71151170575937, -3.35148255039153,$
    -3.01257883277328, -2.67136426606314, -2.3036287214954,$
    -1.87754943767857, -1.32687203939232, -7.72840966450717e-08,$
    1.32687203939232, 1.87754943767857, 2.3036287214954, 2.67136426606314,$
    3.01257883277328, 3.35148255039153, 3.71151170575936, 4.11722526306697,$
    4.59404694329679, 5.16689314460211, 5.85851734698832, 6.68827297583048,$
    7.67162666281268, 8.82031765103987, 10.142824406667, 11.6448185129562,$
    13.3294282264571, 15.1972697734841, 17.2462861630082, 19.4714697638694,$
    21.8645574169981, 24.4137889090064, 27.1038148776609, 29.9158266703621,$
    32.8279553674348, 35.8159497801506, 38.8540980954293, 41.916313892128,$
    44.9772754602266, 48.0134966005524, 51.0042204168578, 53.9320608459731,$
    56.7833601290163, 59.5482728298363, 62.2206194320588, 64.7975706790533,$
    67.2792279882741, 69.6681589022773, 71.9689341802758, 74.1876988925387,$
    76.331796630125, 78.4094552099168, 80.4295344892687, 82.4013318763434,$
    84.3344382773342, 86.2386359278657, 88.123829239849, 90. ]

;mlat = [50.,60.,70.,80.]

;
; FUNCTION EpotVal,lat,mlt,INSIDE=inside,OUTSIDE=outside,FILL_VALUE=fill,BPOT=bpot
;
  nlat = n_elements(mlat)
;
  nmlt = 25
  mlt = fltarr(nmlt)
  for i=0,nmlt-1 do mlt[i] = float(i) ; hourly for one day

; mlt = [0.,6.,12.]
; nmlt = n_elements(mlt)
;
  lon = fltarr(nmlt)
  for i=0,nmlt-1 do lon[i] = mlt[i]*15. ; for plotting
  epot = fltarr(nmlt,nlat)
  fac = fltarr(nmlt,nlat)
;
  tilt = 0.
  swvel = 450.
  swden = 9.
; 
; Hardwire by,bz:
  by = [0.] & bz = [-5.] & bdel = [1.]
; by = [ 0., 0., 0., 0., 0. ,0. ,0. ,0. ,0.]
; bz = [-5.,-4.,-3.,-2.,-1. ,2. ,3. ,4. ,5.]
  nsolar = n_elements(by)
;
; Loop over bmin,bmax for bx,by:
;
; bmin = -10. & bmax = 10. & bdel = 1.
; nsolar = fix((bmax-bmin)/bdel)
; by = fltarr(nsolar) & bz = fltarr(nsolar)
; by[*] = 0.
; for i=0,nsolar-1 do bz[i] = bmin+i*bdel

  if (nsolar ne n_elements(bz)) then begin
    print,'>>> must provide same number of by and bz values:',$
      ' nby=',n_elements(by),' ngz=',n_elements(bz)
    return
  endif
;
; Loop over by,bz conditions:
;
  for isolar = 0,nsolar-1 do begin
;
; Cannot call epotval with bx==0==by
  if (by[isolar] eq 0. and bz[isolar] eq 0.) then begin
    print,'Skipping because by and bz are both zero.'
    continue
  endif 

  setmodel,by[isolar],bz[isolar],tilt,swvel,swden,/yz,$
    path='/oreo/d/foster/weimer05/'

  for j=0,nlat-1 do begin
    for i=0,nmlt-1 do begin
;
; Call EpotVal (w05sc.pro):
      epotv = EpotVal(mlat[j],mlt[i],inside=inside,outside=outside)
      epot[i,j] = epotv

;     if (inside[0] gt 0) then $
;     print,format=$
;       "('test_w05sc: mlat=',f8.3,' mlt=',f8.3,' inside=',i3,' epot=',e12.4)",$
;       mlat[j],mlt[i],inside[0],epot[i,j]
;
; Get field-aligned current from function MPFAC (w05scDerived.pro)
; Note SetMPfac was called from setmodel.
;
      facv = MPFAC(mlat[j],mlt[i],inside=inside,outside=outside)
      fac[i,j] = facv

;     if (inside[0] gt 0) then $
;     print,format=$
;       "('test_w05sc: mlat=',f8.3,' mlt=',f8.3,' inside=',i3,' fac=',e12.4)",$
;       mlat[j],mlt[i],inside[0],fac[i,j]

    endfor ; i=0,nmlt-1
;   if (min(epot[*,j]) ne 0. and max(epot[*,j]) ne 0.) then $
;   print,format="('j=',i3,' mlat=',f8.3,' inside=',i3,' epot min,max=',2e12.4)",$
;     j,mlat[j],inside[0],min(epot[*,j]),max(epot[*,j])
  endfor ; j=0,nlat-1
;
; Contour on stereographic north hemisphere:
;
  if (make_plots eq 0) then continue ; skip plotting
; plot_var = 'EPOT'
  plot_var = 'FAC'
; print,'Plotting ',plot_var
;
  if (keyword_set(ps)) then begin
    set_plot,'PS'
  endif else begin
    set_plot,'X'
    device,retain=2 ; for backing store
  endelse
;
; Center of projection:
  projcen = [90.,0.] ; lat,lon of proj center
  map_stperim = 50.  ; perimeter lat
;
; Position in normalized coords:
;
  xll = .15 & yll = .15 & xur = .85 & yur = .85
  mappos = [xll,yll,xur,yur]
; print,'mappos (xll,yll,xur,yur) =',mappos

  map_set,projcen[0],projcen[1],/stereo,position=mappos,/noborder,$
    limit=[map_stperim,-180.,90.,180.]

  if (plot_var eq 'EPOT') then begin
    fmin = min(epot) & fmax = max(epot)
    print,'epot fmin,max=',fmin,',',fmax
  endif else begin
    fmin = min(fac) & fmax = max(fac)
    print,format="('by=',f8.3,' bz=',f8.3,' fac fmin,max=',e12.4,',',e12.4)",$
      by,bz,fmin,fmax
  endelse
;
; Hardwire contour levels and interval:
;
; cmin = -60. & cmax = 60. & cint = 10. ; hardwire min,max,int
; cmin = fix(fmin) & cmax = fix(fmax)   ; use fmin,fmax
  cmin = fmin & cmax = fmax
  cint = (cmax-cmin)/10.                ; make 10 levels
; cint = 5.                             ; hardwire cint
  nlevels = fix((cmax-cmin)/cint)
  clevels = fltarr(nlevels)
  for i=0,nlevels-1 do clevels[i] = cmin+i*cint
; print,'nlevels=',nlevels,' clevels=' & print,clevels
;
; Negative contours are dashed, positive solid,
; and 0 is long dashes or dotted.
;
  clinestyle = intarr(nlevels)
  clinestyle[*] = 0 ; default solid
  for i=0,nlevels-1 do begin
    if clevels[i] lt 0. then clinestyle[i] = 2 ; dashed
    if clevels[i] eq 0. then clinestyle[i] = 1 ; dotted
;   if clevels[i] eq 0. then clinestyle[i] = 5 ; long dashes
  endfor

  if (plot_var eq 'EPOT') then begin
    contour,epot,lon,mlat,/follow,xstyle=1,ystyle=1,$
      position=mappos,/noerase,/overplot,$
      levels=clevels,c_linestyle=clinestyle
  endif else begin
    contour,fac,lon,mlat,/follow,xstyle=1,ystyle=1,$
      position=mappos,/noerase,/overplot,$
      levels=clevels,c_linestyle=clinestyle
  endelse

  continents = 0
  grid = 0
  if continents eq 1 then map_continents
  if grid eq 1 then map_grid,/label

  labpol_chsize = 1.0
  mag = 1
  hem = 'N'
; fixpolar = 'LON'
  fixpolar = 'SLT'
  labpol,mappos,hem,labpol_chsize,fixpolar,mag=mag

  chsize = 2.0
  if (plot_var eq 'EPOT') then begin
    title = 'Weimer 2005 Electric Potential (kV)'
  endif else begin
    title = 'Weimer 2005 Field-Aligned Current'
  endelse
  xyouts,xll+.5*(xur-xll),yur+.10,title,align=.5,/norm,$
    charsize=chsize

  chsize = 1.2
  title = 'BZ='+string(format="(f8.3)",bz[isolar])+$
         ' BY='+string(format="(f8.3)",by[isolar])+$
         ' SWVEL='+string(format="(f8.3)",swvel)+$
         ' SWDEN='+string(format="(f8.3)",swden)
  xyouts,xll+.5*(xur-xll),yur+.05,title,align=.5,/norm,$
    charsize=chsize

  chsize = 1.2
  minmaxlab = 'MIN,MAX = '+string(format="(g10.4)",fmin)+','+$
                           string(format="(g10.4)",fmax)+'  '+$
             'INTERVAL = '+string(format="(g10.4)",cint)
  xyouts,xll+.5*(xur-xll),yll-.12,minmaxlab,align=.5,/norm,$
    charsize=chsize

  endfor ; isolar=0,nsolar-1

  if (keyword_set(ps)) then device,/close_file

end