load "$NCARG_ROOT/ncl621/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/ncl621/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/ncl621/lib/ncarg/nclscripts/contrib/cd_string.ncl"
load "$NCARG_ROOT/ncl621/lib/ncarg/nclscripts/csm/contributed.ncl"
load "/al11/andrea/research/Elliptical/EllipticalFit.ncl"
;For the dates set below, this code finds the 30km hght contour at 10-hPa and feeds the lat/lon points of the contour
;to a python scripts to calculate the best fit ellipse of that 30 km contour line. (Needed python for imaginary number work)
;This script then reads in the output from that python script (which is in a polar sterographic coord system)
;and converts it into a plotable form to get the lat&lon center, the axes lengths and the angle of rotation of the best fit ellipse 
;This data is then written to the text file EMetricsYYYY1YYYY2
;For 1996: May 1: 122; Jun 1: 153; July 1: 183; Aug 1: 214; Sept 1: 245;  Oct 1: 275; Nov 1: 306; Dec 1:336

begin

lev2plot = 10
lowerlev = 1000
year1 = 1996
year1@calendar = "366_day" ;change to "365_day" for non leap years
mm1 = 11
dd1 = 13

nens = 11
nff = 32

print("y: "+year1+" "+mm1+" "+dd1)

datestr = "1996111300"  ;systemfunc("date -u '+%y%m%d'")

;day1 = day_of_year(year1,mm1,dd1)
date1 = cd_inv_calendar( tointeger(year1), tointeger(mm1), tointeger(dd1), 00, 00, 00,  "hours since 1800-01-01 00:00", 0 )
;date2 = cd_inv_calendar( year1, mm2, dd2, 18, 00, 00,  "hours since 1800-01-01 00:00", 0 )


do ff=0,nff                     ; loop though fcast hrs --> there are ff forecast times

fhr= ff*6                      ; forecast times * 6 gives us forecast hour
fff = sprinti("%0.3i",fhr)     ; string of forecast hour for titles/filenames
if (fhr.lt.100) then           ; check forecast hr to set string correctly for ff time
 fff = sprinti("%0.2i",fhr)    ; gfs reforecast ens data has forecast hours names e.g., 06 09 12 60 108
end if

do p=0,nens-1                   ; loop through ens members --> GEFS reforecast data has 10 ens members and 1 control
pxx = sprinti("%0.2i",p+1)       ; string for ens member number

if (p.eq.0) then               ; when p is 0 - use control run and set some inital stuff up for plotting
 if (ff.eq.0)then              ; first time through do this stuff

  ;grab the analysis data first

  g_file = addfile("/langlab_rit/andrea/Daledata/"+datestr+".p01/gfs_fhr00.p01.grib2","r")  ; Dale data is in langlab_rit
  ; print(getfilevarnames(g_file)) ;;; shows content of g_file

  garray = g_file->HGT_P0_L100_GGA0({lev2plot*100},{0:90},:)   ; var names checked against th print statments above
  uarray = g_file->UGRD_P0_L100_GGA0({lev2plot*100},{65},:)
  tarray = g_file->TMP_P0_L100_GGA0({lev2plot*100},{75:90},:)
  printVarSummary(garray)
  lat = garray&lat_0
  lon = garray&lon_0 
  ylat = dimsizes(lat)
  xlon = dimsizes(lon)

  gfcast = new((/nff+1,nens+1,ylat,xlon/),float,-999)  ; indexes: forecast hr, member, lat, lon
  gfcast(ff,p,:,:)=garray
  ufcast = new((/nff+1,nens+1/),float,-999)
  ufcast(ff,p) = dim_avg(uarray)  ; zonal mean winds at 60
  tfcast = new((/nff+1,nens+1/),float,-999)
  tfcast(ff,p) = dim_avg_wgt(dim_avg(tarray),cos((3.1415/180)*tarray&lat_0),0) ; find the weighted zonal mean avg polar cap t

 end if ; first time through (member 0, ff 0)

 if (ff.ne.0) then  ; 
   print("/langlab_rit/andrea/Daledata/"+datestr+".p"+pxx+"/gfs_fhr"+fff+".p"+pxx+".grib2")
   g_file = addfile("/langlab_rit/andrea/Daledata/"+datestr+".p"+pxx+"/gfs_fhr"+fff+".p"+pxx+".grib2","r")  
   gfcast(ff,p,:,:) = g_file->HGT_P0_L100_GGA0({lev2plot*100},{0:90},:)
   ufcast(ff,p) = dim_avg(g_file->UGRD_P0_L100_GGA0({lev2plot*100},{60},:))
   tfcast(ff,p) = dim_avg_wgt(dim_avg(g_file->TMP_P0_L100_GGA0({lev2plot*100},{75:90},:)),cos((3.1415/180)*tarray&lat_0),0)
      ; find the weighted zonal mean avg polar cap t
 end if
end if ; when p=0
 
if (p.ne.0) then
   g_file = addfile("/langlab_rit/andrea/Daledata/"+datestr+".p"+pxx+"/gfs_fhr"+fff+".p"+pxx+".grib2","r")  
   gfcast(ff,p,:,:) = g_file->HGT_P0_L100_GGA0({lev2plot*100},{0:90},:)
   ufcast(ff,p) = dim_avg(g_file->UGRD_P0_L100_GGA0({lev2plot*100},{60},:))
   tfcast(ff,p) = dim_avg_wgt(dim_avg(g_file->TMP_P0_L100_GGA0({lev2plot*100},{75:90},:)),cos((3.1415/180)*tarray&lat_0),0)
    ; find the weighted zonal mean avg polar cap t
end if

end do ; members
end do ;fcast hrs


printVarSummary(gfcast)
times = date1 ;garray&time
printVarSummary(times)
ntimes = dimsizes(times) 
figstart = 0
print(ntimes)


nmtrc = (nff+1) * (nens+1)
metrics = new(nmtrc,string)

;;;;;;;;;;;Convert times to readable dates;;;;;;;;;;;
 ; Array to hold month abbreviations. Don't store anything in index
  ; '0' (i.e. let index 1=Jan, 2=Feb, ..., index 12=Dec).
  ;
month_abbr = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \
                    "Oct","Nov","Dec"/)
 
 ; Convert to UTC time.

utc_date = cd_calendar(times, 0)
 ; Store return information into more meaningful variables.
 ;
   yy   = tointeger(utc_date(:,0))    ; Convert to integer for
   mm  = tointeger(utc_date(:,1))    ; use sprinti 
   dd    = tointeger(utc_date(:,2))
   hh   = tointeger(utc_date(:,3))

;dayn = day_of_year(yy,mm,dd)

   yystr = sprinti("%0.4i", yy)
   mmstr = month_abbr(mm)
   ddstr = sprinti("%0.2i", dd)
   hhstr = sprinti("%0.2i", hh)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;Begin looping over times;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
do n=0,nff                     ;forecast and analysis times
fignum = figstart + n
print(n)

fhr= n*6 ; 
fff = sprinti("%0.3i",fhr)     ; string of forecast hour for titles/filenames
if (fhr.lt.100) then            ; check forecast hr to set string correctly for ff time
 fff = sprinti("%0.2i",fhr)    ; gfs reforecast ens data has forecast hours names e.g., 06 09 12 60 108
end if

datestr := "F"+fff+": Forecast Initialized: "+hhstr+" UTC "+ddstr+" "+mmstr+" "+yystr

;print(dayn(n))
;doyn = ind(dayind.eq.dayn(n))

do mbr=0,nens

cntr  = new(nens+2,graphic)      

gmean := dim_avg_n_Wrap(gfcast(n,:,:,:),0)

if(mbr.eq.0) then ;--------------start with ens mean plot--------------
;************************************************
; create plot
;************************************************
  wks_type = "png"  
  wks_type@wkWidth = 800
  wks_type@wkHeight = 600 
  filename = "GEFSellipse"+n
  ;filename = "Test"
  wks1 = gsn_open_wks("png","tmp")              ; open a png file
  wks = gsn_open_wks(wks_type ,filename)  
  gsn_define_colormap(wks, "GMT_polar")   ; choose colormap
  res = True 
  res@gsnMaximize = True ;make that map bigger
  res@gsnDraw              = False                ; Don't draw plot 
  res@gsnFrame             = False                ; Don't advance frame


  res@vpHeightF 	 = 0.9               ; change aspect ratio of plot
  res@vpWidthF    = 0.9

 res@gsnRightString = ""       ; turn off right string
  res@gsnLeftString = ""        ; turn off left string

  res@cnFillOn             = False               ; turns on the color
  res@mpFillOn             = False              ; turns off continent gray
  res@cnLinesOn            = True              ; turn off contour lines

  cnres = res    ;make contour specific resource list
  plyres = res   ;make polyline specific resource list
  fres = res

  fres@tiMainString         = lev2plot+"-hPa 30km contour and Best-Fit Ellipse" ; plot title
  fres@tiMainFontHeightF  = 0.021
  fres@gsnCenterString = datestr    ; center title
  fres@gsnCenterStringFontHeightF = 0.019
  fres@cnFillOn          = True              ; color plot desired
  fres@cnLineLabelsOn    = False             ; turn off contour line labels
  fres@cnLineLabelFontHeightF = 0.005
  fres@cnLinesOn         = False             ; turn off contour lines
  fres@gsnSpreadColors     = True            ; use full range of color map
;  fres@cnLevelSelectionMode = "ExplicitLevels"
  fres@cnLevelSelectionMode = "ManualLevels"
  fres@cnInfoLabelOn = False
  fres@cnLevelSpacingF = .25     ; spacing
  fres@cnMinLevelValF = -3.5    ; min value
  fres@cnMaxLevelValF = 3.5    ; max value
;  fres@cnLevels = (/-2,-1.5,-1.,-0.5,0.5,1,1.5,2/)
;  fres@cnMonoLineColor = False
  fres@cnLineThicknessF = 1.5
;  fres@cnLineColors = (/"blue","blue","blue","blue","red","red","red","red"/)
  fres@gsnLeftString = ""
  fres@gsnRightString = ""
; fres@gsnContourNegLineDashPattern = 11
;  fres@gsnContourZeroLineThicknessF = 0.0            ; no zero line
  fres@mpGeophysicalLineColor = "gray80"
  fres@mpGeophysicalLineThicknessF = 2.
  fres@mpProjection      = "Orthographic"       ; choose projection
  fres@mpPerimOn         = False             ; turn off box around plot
  fres@mpFillOn          = False
;  fres@mpMonoFillColor = True
;  fres@mpFillColor = "gray40"
  fres@mpOutlineOn            = True
  fres@mpOutlineDrawOrder     = "PostDraw"
  fres@mpLimitMode            = "LatLon"
  fres@mpMinLatF              =   10.
  fres@mpMaxLatF              =   90.
  fres@mpCenterLonF      = 0.             ; choose center lon
  fres@mpCenterLatF      = 90.               ; choose center lat
  fres@mpGridAndLimbOn   = True              ; turn on lat/lon lines
  fres@mpGridSpacingF    = 10.
  fres@mpGridLineColor       = "gray90"


print("fills")
fills = gsn_csm_map(wks,fres)  ; create the contour plot


 
  cnres@cnLevelSelectionMode = "ExplicitLevels"                       ; set manual/explicit contour levels?
  cnres@cnLevels = (/28500., 29000.,29500.,30000.,30500.,31000./)    ;;;;;;; pick correct contour lines for the level
;  cnres@cnLevels = (/21500.,22000.,22500.,23000.,23500.,24000./)
  cnres@cnMonoLineThickness = True
  cnres@gsLineThicknessesF = (/2.,2.,2., 4., 2., 2./)
  cnres@gsLineThicknessF = 2.
  cnres@cnLineThicknessF = 2.
  cnres@cnLineColor       = "blue" 
  cnres@cnLineLabelFontHeightF = 0.01
  cnres@cnFillDrawOrder      = "Predraw"
  cnres@cnInfoLabelString = ""
  cnres@cnLineLabelBackgroundColor = "transparent" ;the line break at the contour label
  cnres@cnInfoLabelOn = False
;  res@lbLabelBarOn = True 		;Turn on Labelbars
;  res@lbOrientation = "Vertical"	;Label: Veritcal Orientation
;  res@lbLabelPosition = "Right"		;Label: Right Side
;  res@lbBoxMinorExtentF    = .15	;Set Box small length
;  res@lbBottomMarginF = .2		;Set distance btwn edge and labelbar
;  res@lbTopMarginF = .2
;  res@lbLabelFontHeightF = .012		;Set label font height
 cnres@cnLineLabelBackgroundColor = "transparent"

  cnres2=cnres
  cnres2@mpProjection      = "Orthographic"       ; choose projection
  cnres2@mpPerimOn         = False             ; turn off box around plot
  cnres2@mpFillOn          = False
;  cnres2@mpMonoFillColor = True
;  cnres2@mpFillColor = "gray40"
  cnres2@mpOutlineOn            = True
  cnres2@mpOutlineDrawOrder     = "PostDraw"
  cnres2@mpLimitMode            = "LatLon"
  cnres2@mpMinLatF              =   10.
  cnres2@mpMaxLatF              =   90.
  cnres2@mpCenterLonF      = 0.             ; choose center lon
  cnres2@mpCenterLatF      = 90.               ; choose center lat
  cnres2@mpGridAndLimbOn   = True              ; turn on lat/lon lines
  cnres2@mpGridSpacingF    = 10.
  cnres2@mpGridLineColor       = "gray90"


cntr(mbr) = gsn_csm_contour(wks1,gmean,cnres2)    ; create the plot
overlay(fills,cntr(mbr))
print("plotted cnres")
  
; info about contourline 

 isoline := get_isolines(cntr(mbr),30000.)    ;30000 for 10hpa, 23000 for 30hpa, 20000 move for 50
;  cntr(mbr) = gsn_csm_contour(wks,g,cnres)    ; create the plot
; overlay(fills,cntr(mbr))

 plyres@cnLevelSelectionMode = "ManualLevels"
 plyres@gsLineColor       = "mediumblue"                                 ; change to blue
 plyres@gsLineThicknessF = 8.
  count = 0

  segments := isoline@segment_count
  if (segments .gt. 2) then   ;if there are more than 2 its probably an error
    segments := 2               ;so only do 2
  end if                      ;IF USED FOR LEVELS .GT. 50 THIS PROBABLY NEEDS TO BE CHANGED
;  do i = 0,segments -1  
      error = 0   ;reset error to 0 every time
      plyres@gsLineColor       = "mediumblue"                                 ; change to blue every time
;      b := isoline@start_point(i)
;      e := b + isoline@n_points(i) - 1
      ylat := isoline(0,:)             ; b:e)  <--------------- let try reading all points in 
      xlon := isoline(1,:)             ; b:e)  <--------------- and fitting ellipse to all points
   
 if (dimsizes(ylat) .le. 50) then    ;it its an error set error = 1
    error = 1
 end if


;;;;;;;remove comments of lines with i for doing segemnt things
; if (i .eq. 0 .and. error .eq. 0) then    ;there's only 1 vortex or its the first one
      plyplot1 = gsn_add_polyline(wks,fills,xlon,ylat,plyres)
; end if
      
; if (i .eq. 1 .and. error .eq. 0) then ;it really is a split vortex
;      plyplot2 = gsn_add_polyline(wks,fills,xlon,ylat,plyres)
; end if
 
; count = count + isoline@n_points(i)
; print(isoline@level + " has " + count + " total points in " + isoline@segment_count + " segments" )
max_npts = max(isoline@n_points)

lats := ylat*(3.1415/180)
lons := xlon*(3.1415/180)

;convert lat and lon to new coords with NP at origin and 0 lon as y
; Following Waugh (1997)
;;;;;;;Write contour lat/lon to file;;;;;;;;;;;;;;;
x := (cos(lons) * cos(lats))/(1 + sin(lats))
asciiwrite("x.txt",x)
y := (sin(lons) * cos(lats))/(1 + sin(lats))
asciiwrite("y.txt",y)

;;;;;;;;Run Python code to find Ellipse;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;uses x.txt and y.txt to find ellipse then makes xx.txt and yy.txt;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
cmd = "python fitEllipse2.py"
system(cmd)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
print("Ran Python")
xx := asciiread("xx.txt",-1,"float")
yy := asciiread("yy.txt",-1,"float")

;convert polar sterographic cartesian coords to lat/lon
xxlon := atan(yy/xx)
xxlon := where(xx.lt.0, where(yy.gt.0, atan(yy/xx) +3.1415 ,atan(yy/xx)-3.14159), atan(yy/xx))
yysinxxlon := yy / sin(xxlon)
pi4 = 3.1415 / 4
yylat := -2 * ( atan(yysinxxlon) - pi4)  
lats := yylat*(180/3.1415)
lons := xxlon*(180/3.1415)

do l=1,dimsizes(lats)-1
 if(abs(lats(l)-lats(l-1)).gt.1.5)then
  lats(l) = lats(l-1)
 end if
end do

;print(lats)
;print(lons)


;read in major and minor axis lengths in cart. coords, center of ellipse, and rotation of major axis in radians
vals := asciiread("axescenterphi.txt", -1,"float")
print(vals)
aax := vals(0)
bax := vals(1)

centerx := vals(2)
centery := vals(3)
phi := vals(4)

;convert cartisian to lat/lon coords
; if x = ( cos(lon)cos(lat) ) / (1 + sin(lat)) and y = ( sin(lon)cos(lat) ) / (1 + sin(lat))
; then lon = atan(y/x) and lat = 2[atan( y/sin(lon) ) - pi/4]  
;      via tangent half angle formula (see http://en.wikipedia.org/wiki/Tangent_half-angle_formula )
cenlon := atan(centery / centerx) 
cenlon := where(centerx.lt.0, where(centery.gt.0, atan(centery/centerx) +3.1415 ,atan(centery/centerx)-3.14159), atan(centery/centerx))
cenlondeg := cenlon*180/3.1415
ysinlon := centery / sin(cenlon)
cenlat := -2 * ( atan(ysinlon) - pi4)  
cenlatdeg := cenlat*(180/3.1415)

print("center of ellipse = "+cenlatdeg+" N, "+cenlondeg+" E.")


;Calc enpoints of the axes of the vortex 
xa = (aax) * cos(phi)
ya = (aax) * sin(phi)
xb = (bax) * sin(phi)
yb = (bax) * cos(phi)

endx = (/ centerx+xa , centerx-xa , centerx+xb , centerx-xb/) 
endy = (/ centery+ya , centery-ya , centery-yb , centery+yb/)
endlon = atan(endy / endx)
endlon = where(endx.lt.0, where(endy.gt.0, atan(endy/endx) +3.1415 ,atan(endy/endx)-3.14159), atan(endy/endx))
endlondeg = endlon*180/3.1415
ysinelon = endy / sin(endlon)
endlat = -2 * ( atan(ysinelon) - pi4)  
endlatdeg = endlat*(180/3.1415)
;print(endlatdeg+" "+endlondeg)

;add lines for the a and b ellipse axes
  resp                  = True                      ; polyline mods desired
  resp@gsLineColor      = "firebrick"                     ; color of lines
  resp@gsLineThicknessF = 2.0                       ; thickness of lines

; create array of dummy graphic variables. This is required, b/c each line
; must be associated with a unique dummy variable.

; find distance and points on the great circles that represent the major and minor ellipse axes
  a1gcircle := gc_latlon(endlatdeg(0),endlondeg(0),cenlatdeg,cenlondeg,100,4) 
  a2gcircle := gc_latlon(cenlatdeg,cenlondeg,endlatdeg(1),endlondeg(1),100,4) 
  b1gcircle := gc_latlon(endlatdeg(2),endlondeg(2),cenlatdeg,cenlondeg,100,4) 
  b2gcircle := gc_latlon(cenlatdeg,cenlondeg,endlatdeg(3),endlondeg(3),100,4) 
;  dum := new(4,graphic)  

; draw each line separately. Each line must contain two points.
 resp@gsLineLabelString= a1gcircle+" km"                    ; adds a line label string
;   dum(0)=gsn_add_polyline(wks,fills,a1gcircle@gclon,a1gcircle@gclat,resp)      ;a1 axis
;   dum(1)=gsn_add_polyline(wks,fills,a2gcircle@gclon,a2gcircle@gclat,resp)      ;a1 axis
 resp@gsLineLabelString= b1gcircle+" km"                    ; adds a line label string
;   dum(2)=gsn_add_polyline(wks,fills,b1gcircle@gclon,b1gcircle@gclat,resp)      ;b axis   
;   dum(3)=gsn_add_polyline(wks,fills,b2gcircle@gclon,b2gcircle@gclat,resp)      ;b axis   




plyres@gsLineColor       = "deepskyblue"                                 ; change to blue
ellipseplot = gsn_add_polyline(wks,fills,lons,lats,plyres)
overlay(fills,ellipseplot)

plyres@gsMarkerColor = "deepskyblue"
plyres@gsMarkerIndex = 1
plyres@gsMarkerSizeF = 0.09
ellipseplot2 = gsn_add_polymarker(wks,fills,cenlondeg,cenlatdeg,plyres)
overlay(fills,ellipseplot2)

 

;; make sure we write major axis first and that phi corresponds to angle of major axis
phideg = phi * 180./3.1415  ; if a1 is greater than b1
if (a1gcircle .lt. b1gcircle) then  ; if a1 is less than b1, adjust phi by 90 deg
  phideg = phi * 180./3.1415 - 90
end if

print("Emetrics phi: "+phideg)

newphideg = phideg
if (phideg.lt.-45) then
  newphideg = phideg+180
end if

ratio = a1gcircle/b1gcircle
if (ratio.lt.1) then
ratio = 1/ratio
end if

size=3.14159*a1gcircle*b1gcircle

Ulabel = sprintf("%5.2f",dim_avg_Wrap(ufcast(n,:)))
Tlabel = sprintf("%5.2f",dim_avg_Wrap(tfcast(n,:)))
slabel = sprintf("%5.2e",size)
philabel = sprintf("%5.2f",newphideg)
rlabel = sprintf("%5.2f",ratio)
lonlabel = sprintf("%5.2f",cenlondeg)

;;;;Create legend with Elliptical Diagnostics
   lgres                    = True
   lgres@lgLineColors       = (/"white","white","white","white", "white", "white","white","white"/)
   lgres@lgLineThicknessF   = 1. 
   lgres@lgLabelFontHeightF = .38            ; set the legend label font thickness
   lgres@vpWidthF           = 0.25           ; width of legend (NDC)
   lgres@vpHeightF          = 0.23            ; height of legend (NDC)
   lgres@lgMonoDashIndex    = True   
   lgres@lgPerimColor       = "deepskyblue"       ; draw the box perimeter in orange
  ;lgres@lgPerimThicknessF = 2.0            ; thicken the box perimeter
   lgres@lgPerimFill = 0   ; solid fill
   lgres@lgPerimFillColor = "White"
   labels = (/ "U65: "+Ulabel, "Polar Cap T: "+Tlabel, "Size: "+slabel, "Rotation: "+philabel+"~S~o~N~","Ratio: "+rlabel, "Center Lon: "+lonlabel+"~S~o~N~E ", "Center Lat: "+sprintf("%5.2f",cenlatdeg)+"~S~o~N~N ","~F22~Control Run" /)

; Create the legend.
;   lbid   = gsn_create_legend(wks,8,labels,lgres)         ; create legend   <--- Control Run box off

; Set up resources to attach legend to map.
   amres = True
   amres@amParallelPosF   =  0.4 	         ; positive move legend to the right
   amres@amOrthogonalPosF =  0.4                 ; positive move the legend up
;   annoid1 = gsn_add_annotation(fills,lbid,amres)   ; attach legend to plot

enssize := new(nens+1,float,0)    ; 0 -> ens mean, >0 ens members
ensphi := new(nens+1,float,0)
enslat := new(nens+1,float,0)
enslon := new(nens+1,float,0)
ensratio := new(nens+1,float,0)
ensu := new(nens+1,float,0)
enst := new(nens+1,float,0)
ensseg := new(nens+1,float,0)

;if (isoline@n_points(i).eq.max_npts) then   <-------add back in for segment looping
ensseg(0) = segments
enssize(0) = size
ensphi(0) = newphideg
enslat(0) = cenlatdeg
enslon(0) = cenlondeg
ensratio(0) = ratio
aaa = a1gcircle
bbb = b1gcircle
ensu(0) = dim_avg_Wrap(ufcast(n,:))  ;<----ens mean info
enst(0) = dim_avg_Wrap(tfcast(n,:))  ;<----ens mean info
;end if                                <-------add back in for segment looping

;  end do; end do for segments in member = 0  <-------------------add back in for segments loops

end if ;when mbr=0---------------------------------------------------
;---------------------------------------------------------------------------------------------


;------------------------Plotting for ENS members---------------------------
;---------------------------------------------------------------------------------------------

if (mbr.ne.0) then  ;;;;; this loop is for ens members
g := gfcast(n,mbr-1,:,:)             ; why mbr-1, b/c in the first iteration of the mbr loop we plot the mbr mean, 
                                     ; second itteration starts the plotting of indivd. members

cnres@cnLevels = (/2850., 2900.,2950.,3000.,3050.,3100./)
printVarSummary(g)
  cntr(mbr) = gsn_csm_contour(wks1,g,cnres2)    ; create the plot
  
print("plotted cnres")
  
; info about contourline 

 isoline := get_isolines(cntr(mbr),30000.)    ;30000 for 10hpa, 23000 for 30hpa, 20000 move for 50
  cntr(mbr) = gsn_csm_contour(wks,g,cnres)    ; create the plot
overlay(fills,cntr(mbr))

 plyres@cnLevelSelectionMode = "ManualLevels"
 plyres@gsLineColor       = -1                               ; change to transparent
 plyres@gsLineThicknessF = 8.
  count = 0
  segments := 0
  segments := isoline@segment_count
  if (segments .gt. 2) then   ;if there are more than 2 its probably an error
    segments := 2               ;so only do 2
  end if                      ;IF USED FOR LEVELS .GT. 50 THIS PROBABLY NEEDS TO BE CHANGED
;  do i = 0,segments -1             ; <---- add back in for known splits
      error = 0   ;reset error to 0 every time
      plyres@gsLineColor       = "mediumblue"                                 ; change to blue every time
;     b := isoline@start_point(i)           ; <---- add back in for known splits
;      e := b + isoline@n_points(i) - 1     ; <---- add back in for known splits
      ylat := isoline(0,:)   ; b:e)     ; <---- add back in for known splits
      xlon := isoline(1,:)   ; b:e)     ; <---- add back in for known splits
   
 if (dimsizes(ylat) .le. 50) then    ;it its an error set error = 1
    error = 1
 end if

; if (i .eq. 0 .and. error .eq. 0) then    ;there's only 1 vortex or its the first one  ; <---- add back in for known splits
      plyplot1 = gsn_add_polyline(wks,fills,xlon,ylat,plyres)
; end if                                                   ; <---- add back in for known splits
      
; if (i .eq. 1 .and. error .eq. 0) then ;it really is a split vortex    ; <---- add back in for known splits
;      plyplot2 = gsn_add_polyline(wks,fills,xlon,ylat,plyres)          ; <---- add back in for known splits
; end if                                                                ; <---- add back in for known splits
 
; count = count + isoline@n_points(i)
; print(isoline@level + " has " + count + " total points in " + isoline@segment_count + " segments" )

;  printVarSummary(xlon)
;  print("Avg lon: "+dim_avg(xlon))
;  printVarSummary(ylat)
;  print("Avg lat: "+dim_avg(ylat))
;  print(x)
;  print(y)


lats := ylat*(3.1415/180)
lons := xlon*(3.1415/180)

;convert lat and lon to new coords with NP at origin and 0 lon as y
; Following Waugh (1997)
;;;;;;;Write contour lat/lon to file;;;;;;;;;;;;;;;
x := (cos(lons) * cos(lats))/(1 + sin(lats))
asciiwrite("x.txt",x)
y := (sin(lons) * cos(lats))/(1 + sin(lats))
asciiwrite("y.txt",y)

;;;;;;;;Run Python code to find Ellipse;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;uses x.txt and y.txt to find ellipse then makes xx.txt and yy.txt;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
cmd = "python fitEllipse2.py"
system(cmd)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 
print("Ran Python")
xx := asciiread("xx.txt",-1,"float")
yy := asciiread("yy.txt",-1,"float")

;convert polar sterographic cartesian coords to lat/lon
xxlon := atan(yy/xx)
xxlon := where(xx.lt.0, where(yy.gt.0, atan(yy/xx) +3.1415 ,atan(yy/xx)-3.14159), atan(yy/xx))
sinxxlon = sin(xxlon)
sinxxlon  = where(sinxxlon.eq.0, sinxxlon@_FillValue,sinxxlon)   ; don't divide by 0 

 yysinxxlon := yy / sinxxlon

pi4 = 3.1415 / 4
yylat := -2 * ( atan(yysinxxlon) - pi4)  
lats := yylat*(180/3.1415)
lons := xxlon*(180/3.1415)


do l=1,dimsizes(lats)-1
 if(abs(lats(l)-lats(l-1)).gt.1.5)then
  lats(l) = lats(l-1)
 end if
end do


;read in major and minor axis lengths in cart. coords, center of ellipse, and rotation of major axis in radians
vals := asciiread("axescenterphi.txt", -1,"float")
print(vals)
aax := vals(0)
bax := vals(1)

centerx := vals(2)
centery := vals(3)
phi := vals(4)

;convert cartisian to lat/lon coords
; if x = ( cos(lon)cos(lat) ) / (1 + sin(lat)) and y = ( sin(lon)cos(lat) ) / (1 + sin(lat))
; then lon = atan(y/x) and lat = 2[atan( y/sin(lon) ) - pi/4]  
;      via tangent half angle formula (see http://en.wikipedia.org/wiki/Tangent_half-angle_formula )
cenlon := atan(centery / centerx) 
cenlon := where(centerx.lt.0, where(centery.gt.0, atan(centery/centerx) +3.1415 ,atan(centery/centerx)-3.14159), atan(centery/centerx))
cenlondeg := cenlon*180/3.1415
ysinlon := centery / sin(cenlon)
cenlat := -2 * ( atan(ysinlon) - pi4)  
cenlatdeg := cenlat*(180/3.1415)

print("center of ellipse = "+cenlatdeg+" N, "+cenlondeg+" E.")


;Calc enpoints of the axes of the vortex 
xa = (aax) * cos(phi)
ya = (aax) * sin(phi)
xb = (bax) * sin(phi)
yb = (bax) * cos(phi)

endx = (/ centerx+xa , centerx-xa , centerx+xb , centerx-xb/) 
endy = (/ centery+ya , centery-ya , centery-yb , centery+yb/)
endlon = atan(endy / endx)
endlon = where(endx.lt.0, where(endy.gt.0, atan(endy/endx) +3.1415 ,atan(endy/endx)-3.14159), atan(endy/endx))
endlondeg = endlon*180/3.1415
ysinelon = endy / sin(endlon)
endlat = -2 * ( atan(ysinelon) - pi4)  
endlatdeg = endlat*(180/3.1415)
;print(endlatdeg+" "+endlondeg)

;add lines for the a and b ellipse axes
  resp                  = True                      ; polyline mods desired
  resp@gsLineColor      = "firebrick"                     ; color of lines
  resp@gsLineThicknessF = 2.0                       ; thickness of lines

; create array of dummy graphic variables. This is required, b/c each line
; must be associated with a unique dummy variable.

; find distance and points on the great circles that represent the major and minor ellipse axes
  a1gcircle := gc_latlon(endlatdeg(0),endlondeg(0),cenlatdeg,cenlondeg,100,4) 
  a2gcircle := gc_latlon(cenlatdeg,cenlondeg,endlatdeg(1),endlondeg(1),100,4) 
  b1gcircle := gc_latlon(endlatdeg(2),endlondeg(2),cenlatdeg,cenlondeg,100,4) 
  b2gcircle := gc_latlon(cenlatdeg,cenlondeg,endlatdeg(3),endlondeg(3),100,4) 
;  dum := new(4,graphic)  

; draw each line separately. Each line must contain two points.
; resp@gsLineLabelString= a1gcircle+" km"                    ; adds a line label string
;   dum(0)=gsn_add_polyline(wks,fills,a1gcircle@gclon,a1gcircle@gclat,resp)      ;a1 axis
;   dum(1)=gsn_add_polyline(wks,fills,a2gcircle@gclon,a2gcircle@gclat,resp)      ;a1 axis
; resp@gsLineLabelString= b1gcircle+" km"                    ; adds a line label string
;   dum(2)=gsn_add_polyline(wks,fills,b1gcircle@gclon,b1gcircle@gclat,resp)      ;b axis   
;   dum(3)=gsn_add_polyline(wks,fills,b2gcircle@gclon,b2gcircle@gclat,resp)      ;b axis   


;if (i .eq. 0 .and. error .eq. 0) then              ; <---- add back in for known splits
plyres@gsLineThicknessF = 2.0                       ; thickness of lines
plyres@gsLineColor       = "slategray"                                 ; change to blue
 if (segments.eq.2) then
  plyres@gsLineColor       = "lavenderblush3"                                 ; if split change line color
 end if
ellipseplot := gsn_add_polyline(wks,cntr(mbr),lons,lats,plyres)
overlay(fills,ellipseplot)



;plyres@gsMarkerColor = "lightslategray"
;plyres@gsMarkerIndex = 4
;plyres@gsMarkerSizeF = 0.01
;ellipseplot2 := gsn_add_polymarker(wks,cntr(mbr),cenlondeg,cenlatdeg,plyres)
;end if                                           ; <---- add back in for known splits
 
;if (i .eq. 1 .and. error .eq. 0) then            ; <---- add if statement back in for known splits

;plyres@gsLineThicknessF = 2.0                       ; thickness of lines
;plyres@gsLineColor       = "slategray"                                 ; change to blue
; if (segments.eq.2) then
;   plyres@gsLineColor       = "lavenderblush3"                                 ; if split change line color
; end if
;ellipseplot2 := gsn_add_polyline(wks,cntr(mbr),lons,lats,plyres)
;overlay(fills,ellipseplot2)

;end if



;; make sure we write major axis first and that phi corresponds to angle of major axis
phideg = phi * 180./3.1415  ; if a1 is greater than b1
if (a1gcircle .lt. b1gcircle) then  ; if a1 is less than b1, adjust phi by 90 deg
  phideg = phi * 180./3.1415 - 90
end if

print("Emetrics phi: "+phideg)

newphideg = phideg
if (phideg.lt.-45) then
  newphideg = phideg+180
end if

ratio = a1gcircle/b1gcircle
if (ratio.lt.1) then
ratio = 1/ratio
end if

size=3.14159*a1gcircle*b1gcircle
max_npts = max(isoline@n_points)

print("n: "+n+",  mbr: "+mbr)
;if (isoline@n_points(i).eq.max_npts) then         ; <---- add back in for known splits
ensseg(mbr) = segments
enssize(mbr) = size
ensphi(mbr) = newphideg
enslat(mbr) = cenlatdeg
enslon(mbr) = cenlondeg
aaa = a1gcircle
bbb = b1gcircle
ensratio(mbr) = ratio
ensu(mbr) = ufcast(n,mbr)
enst(mbr) = tfcast(n,mbr)
; end if                                        ; <---- add back in for known splits

if (mbr.eq.11) then
 
do w=0,dimsizes(enslat)-1
 print("In loop at L743, w = "+w)
 plyres@gsMarkerColor = "slategray"
; if (ensseg(mbr).gt.1) then
;   plyres@gsMarkerColor       = "lavenderblush3"    ; if split change line color  <-----add for segment 
; end if
 plyres@gsMarkerIndex = 4
 plyres@gsMarkerSizeF = 0.01
 ellipseplot2 := gsn_add_polymarker(wks,cntr(mbr),enslon(w),enslat(w),plyres)
 overlay(fills,ellipseplot2)
end do

;;;;Create legend with Elliptical Diagnostics
   lgres                    = True
   lgres@lgLineColors       = (/"white","white","white","white", "white", "white","white","white"/)
   lgres@lgLineThicknessF   = 1. 
   lgres@lgLabelFontHeightF = .38            ; set the legend label font thickness
   lgres@vpWidthF           = 0.45           ; width of legend (NDC)
   lgres@vpHeightF          = 0.23            ; height of legend (NDC)
   lgres@lgMonoDashIndex    = True   
   lgres@lgPerimColor       = "slategray"       ; draw the box perimeter in orange
  ;lgres@lgPerimThicknessF = 3.0            ; thicken the box perimeter
   lgres@lgPerimFill = 0   ; solid fill
   lgres@lgPerimFillColor = "White"
   labels = (/"U65: "+sprintf("%5.2f",dim_avg(ensu(:)) )+", "+sprintf("%5.2f",min(ensu(:)) )+", "+sprintf("%5.2f",max(ensu(:)) ),"Polar Cap T: "+sprintf("%5.2f",dim_avg(enst(:)) )+", "+sprintf("%5.2f",min(enst(:)) )+", "+sprintf("%5.2f",max(enst(:)) ),"Size: "+sprintf("%5.2e",dim_avg(enssize(:)))+", "+sprintf("%5.2e",min(enssize(:)))+", "+sprintf("%5.2e",max(enssize(:))),"Rotation: "+sprintf("%5.2f",dim_avg(ensphi(:)))+"~S~o~N~,"+sprintf("%5.2f",min(ensphi(:)))+"~S~o~N~,"+sprintf("%5.2f",max(ensphi(:)))+"~S~o~N~", "Ratio: "+sprintf("%5.2f",dim_avg(ensratio(:)))+", "+sprintf("%5.2f",min(ensratio(:)))+", "+sprintf("%5.2f",max(ensratio(:))), "Center Lon: "+sprintf("%5.2f",dim_avg(enslon(:)))+"~S~o~N~E, "+sprintf("%5.2f",min(enslon(:)))+"~S~o~N~E, "+sprintf("%5.2f",max(enslon(:)))+"~S~o~N~E ", "Center Lat: "+sprintf("%5.2f",dim_avg(enslat(:)))+"~S~o~N~N, "+sprintf("%5.2f",min(enslat(:)))+"~S~o~N~N, "+sprintf("%5.2f",max(enslat(:)))+"~S~o~N~N ","~F22~Ens. Mean/Min/Max" /)

; Create the legend.
   lbid   = gsn_create_legend(wks,8,labels,lgres)         ; create legend

; Set up resources to attach legend to map.
   amres = True
   amres@amParallelPosF   =  0.3 	         ; positive move legend to the right
   amres@amOrthogonalPosF =  -0.37                 ; positive move the legend up
   annoid1 = gsn_add_annotation(fills,lbid,amres)   ; attach legend to plot
end if ;mbr=11


; end do  ;ending segment do      ; <---- add back in for known splits
end if ;-------------------------------------------------------------


nmbr = n*nens + mbr     

;;;;write the Elliptical diagnostics to a file;;;;;;;
metrics(nmbr) =  hhstr+" "+ddstr+" "+mmstr+" "+yystr+" "+fff+" "+mbr+" "+enslon(mbr)+" "+enslat(mbr)+" "+aaa+" "+bbb+" "+ensphi(mbr)+" "+ufcast(n,mbr)+" "+tfcast(n,mbr)+" "+ensseg(mbr)




end do; loop through members



draw(fills)

 pres = True
 maximize_output(wks,pres)
 delete(wks1)

end do
;;;;;;;;;;;;;;;;;;;;;;;;;end the time do loop;;;;;;;;;

asciiwrite("GEFSmetrics"+year1+mm1+dd1+".txt", metrics)

end