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 = 30
lowerlev = 1000
year1 = systemfunc("date -u '+%Y'")
year1@calendar = "365_day" ;change to "365_day" for non leap years
mm1 = systemfunc("date -u '+%m'")
dd1 = systemfunc("date -u '+%d'")
mm2 = 11
dd2 = 30
print("y: "+year1+" "+mm1+" "+dd1)

datestr = 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,64 ;  run though fcast hrs

fhr= ff*6 ; 
fff = sprinti("%0.3i",fhr)

do p=0,20  ; run through ens members
pxx = sprinti("%0.2i",p)

if (p.eq.0) then ; use control run
 if (ff.eq.0)then  ; first time through
  ;grab the analysis data first
  g_file = addfile("/cas2/unidata/GRIB/gefs/GFS_"+datestr+"_c00_00_000.grb2","r")  
  print(getfilevarnames(g_file))

  garray = g_file->HGT_P1_L100_GLL0({lev2plot*100},{0:90},:)
  uarray = g_file->UGRD_P1_L100_GLL0({lev2plot*100},{65},:)
  tarray = g_file->TMP_P1_L100_GLL0({lev2plot*100},{75:90},:)
  printVarSummary(garray)
  lat = garray&lat_0
  lon = garray&lon_0 
  ylat = dimsizes(lat)
  xlon = dimsizes(lon)

  gfcast = new((/65,21,ylat,xlon/),float,-999)  ; indexes: forecast hr, member, lat, lon
  gfcast(ff,p,:,:)=garray
  ufcast = new((/65,21/),float,-999)
  ufcast(ff,p) = dim_avg(uarray)  ; zonal mean winds at 60
  tfcast = new((/65,21/),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  ; 
   g_file = addfile("/cas2/unidata/GRIB/gefs/GFS_"+datestr+"_c"+pxx+"_00_"+fff+".grb2","r")  
   gfcast(ff,p,:,:) = g_file->HGT_P1_L100_GLL0({lev2plot*100},{0:90},:)
   ufcast(ff,p) = dim_avg(g_file->UGRD_P1_L100_GLL0({lev2plot*100},{60},:))
   tfcast(ff,p) = dim_avg_wgt(dim_avg(g_file->TMP_P1_L100_GLL0({lev2plot*100},{75:90},:)),cos((3.1415/180)*tarray&lat_0),0)
 end if
end if ; when p=0
 
if (p.ne.0) then
   g_file = addfile("/cas2/unidata/GRIB/gefs/GFS_"+datestr+"_p"+pxx+"_00_"+fff+".grb2","r")  
   gfcast(ff,p,:,:) = g_file->HGT_P1_L100_GLL0({lev2plot*100},{0:90},:)
   ufcast(ff,p) = dim_avg(g_file->UGRD_P1_L100_GLL0({lev2plot*100},{60},:))
   tfcast(ff,p) = dim_avg_wgt(dim_avg(g_file->TMP_P1_L100_GLL0({lev2plot*100},{75:90},:)),cos((3.1415/180)*tarray&lat_0),0)
end if

end do ; members
end do ;fcast hrs


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



metrics = new(65*21,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,64 ;forecast and analysis times
fignum = figstart + n
print(n)

fhr= n*6 ; 
fff = sprinti("%0.3i",fhr)
datestr := "F"+fff+": Forecast Initialized: "+hhstr+" UTC "+ddstr+" "+mmstr+" "+yystr

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

do mbr=0,20

cntr  = new(21,graphic)      
g := gfcast(n,mbr,:,:)

if(mbr.eq.0) then ;--------------start with control run--------------
;************************************************
; create plot
;************************************************
  wks_type = "png"  
  wks_type@wkWidth = 800
  wks_type@wkHeight = 600 
  filename = "GEFS/GEFSellipse30"+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./)
  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       = "darkblue" 
  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,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       = "mediumblue"                                 ; change to blue
 plyres@gsLineThicknessF = 8.
  count = 0
  do i = 0, isoline@segment_count -1
      b := isoline@start_point(i)
      e := b + isoline@n_points(i) - 1
      ylat := isoline(0,b:e)
      xlon := isoline(1,b:e)
      plyplot = gsn_add_polyline(wks,fills,xlon,ylat,plyres)
      count = count + isoline@n_points(i)
  end do

  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))
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",ufcast(n,mbr))
Tlabel = sprintf("%5.2f",tfcast(n,mbr))
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

; 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(20,float,0)
ensphi := new(20,float,0)
enslat := new(20,float,0)
enslon := new(20,float,0)
ensratio := new(20,float,0)
ensu := new(20,float,0)
enst := new(20,float,0)

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


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

if (mbr.ne.0) then

cnres@cnLevels = (/2850., 2900.,2950.,3000.,3050.,3100./)

  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
  do i = 0, isoline@segment_count -1
      b := isoline@start_point(i)
      e := b + isoline@n_points(i) - 1
      ylat := isoline(0,b:e)
      xlon := isoline(1,b:e)
      plyplot = gsn_add_polyline(wks,fills,xlon,ylat,plyres)
      count = count + isoline@n_points(i)
  end do

  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))
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


;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@gsLineThicknessF = 2.0                       ; thickness of lines
plyres@gsLineColor       = "slategray"                                 ; change to blue
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)

 

;; 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

print("n: "+n+",  mbr: "+mbr)
enssize(mbr-1) = size
ensphi(mbr-1) = newphideg
enslat(mbr-1) = cenlatdeg
enslon(mbr-1) = cenlondeg
ensratio(mbr-1) = ratio
ensu(mbr-1) = ufcast(n,mbr)
enst(mbr-1) = tfcast(n,mbr)

if (mbr.eq.20)then
 
do w=0,dimsizes(enslat)-1
 plyres@gsMarkerColor = "slategray"
 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.25           ; 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(:)) ),"Polar Cap T: "+sprintf("%5.2f",dim_avg(enst(:)) ),"Size: "+sprintf("%5.2e",dim_avg(enssize(:))),"Rotation: "+sprintf("%5.2f",dim_avg(ensphi(:)))+"~S~o~N~", "Ratio: "+sprintf("%5.2f",dim_avg(ensratio(:))), "Center Lon: "+sprintf("%5.2f",dim_avg(enslon(:)))+"~S~o~N~E ", "Center Lat: "+sprintf("%5.2f",dim_avg(enslat(:)))+"~S~o~N~N ","~F22~Ens. Mean" /)

; 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.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
end if ;mbr=20


end if ;-------------------------------------------------------------



nmbr = n*20+mbr

;;;;write the Elliptical diagnostics to a file;;;;;;;
metrics(nmbr) =  hhstr+" "+ddstr+" "+mmstr+" "+yystr+" "+fff+" "+mbr+" "+cenlondeg+" "+cenlatdeg+" "+a1gcircle+" "+b1gcircle+" "+phideg+" "+ufcast(n,mbr)+" "+tfcast(n,mbr)







end do; loop through members


draw(fills)

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

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

asciiwrite("GEFS/GEFSmetrics30"+year1+mm1+dd1+".txt", metrics)

end