;+
; NAME:
;   PLOT_FLEXION
;
;
; AUTHOR:
;   Benoit Carry, ESO / LESIA, Observatoire de Paris
;
;
; PURPOSE:
;
;
; CALLING SEQUENCE:
;
;
; INPUT:
;   MAP - Map to be plotted
;
;
; OUTPUT:
;
;
; KEYWORD PARAMETERS:
;
;
;
; CREATION
;   2007-Nov-20: B.Carry (ESO), based on Goldeberg-Gott routines arXiv 2006
;
;
; HISTORY
;   2007-Jan-11: B.Carry (ESO) Keyword update: use of structure
;
;-
pro plotMapFlexion, lat, lon, projection=projection, flexINFO=flexINFO



;----------------------------------------;
;- Initilization And Input Verification -;
;----------------------------------------;

  if not keyword_set(projection) then begin
    print, 'SPECIFY A PROJECTION'
    stop
  endif


 ;-Mapping Deformation
  if not keyword_set(flexionINFO) then $
     flexINFO={latStart:-80  ,$    ;-Flexion: Latitude lower value  
               latEnd  : 80  ,$    ;-Flexion: Latitude upper value  
               latStep : 20  ,$    ;-Flexion: Latitude step         
               longStart: 45 ,$    ;-Flexion: Longitude lower value 
               longEnd  : 405,$    ;-Flexion: Longitude upper value 
               longStep : 90 ,$    ;-Flexion: Longitude step        
               Radius: 6     }     ;-Flexion: Circle radius


linecolor=0



;--------------------------------------;
;- Projection Deformation Calculation -;
;--------------------------------------;

 ;-Ellipse smoothness
  ntheta=400
  theta=indgen(ntheta)*2*!DPI/float(ntheta)


npts=200

dlon0=fltarr(npts)
dlat0=fltarr(npts)

ux_arr=fltarr(npts)
uy_arr=fltarr(npts)

  for lat=flexINFO.latStart, flexINFO.latEnd, flexINFO.latStep do $
     for lon=flexINFO.longStart,flexINFO.longEnd, flexINFO.longStep do begin

if lat eq 0 then linecolor=255
if lat ne 0 then linecolor=0

 ;-Flexion Computation
  ux=1./cos(lat*!DTOR)
  uy=0
  dlon0[0]=0
  dlat0[0]=0
  dtau=4.*flexINFO.Radius/float(npts)
  for i=1,npts/4-1 do begin
    dlon0[i]=dlon0[i-1]+ux*dtau
    dlat0[i]=dlat0[i-1]+uy*dtau
    ux=ux+2*tan((lat+dlat0[i])*!DTOR)*!DTOR*ux*uy*dtau
    uy=uy-0.5*sin(2*(lat+dlat0[i])*!DTOR)*!DTOR*ux*ux*dtau
  endfor
  ux=-1./cos(lat*!DTOR)
  uy=0
  dlon0[npts/4]=0
  dlat0[npts/4]=0
  for i=npts/4+1,npts/2-1 do begin
    dlon0[i]=dlon0[i-1]+ux*dtau
    dlat0[i]=dlat0[i-1]+uy*dtau
    ux=ux+2*tan((lat+dlat0[i])*!DTOR)*!DTOR*ux*uy*dtau
    uy=uy-0.5*sin(2*(lat+dlat0[i])*!DTOR)*!DTOR*ux*ux*dtau
  endfor
  for i=0,npts/2-1 do begin
    dlon0[i+npts/2]=0
    dlat0[i+npts/2]=-flexINFO.Radius+2.*flexINFO.Radius*(i+0.5)/(npts/2)
  endfor
  deproject,x,y,lat+dlat0,lon+dlon0,projection=projection,forward=1
  x = x mod (2.*!DPI)
            



;------------------;
;- Cross Plotting -;
;------------------;

 ;-Cross Horizontal Axis - Right
  idx=where((lon+dlon0)[0:npts/4-1] ge flexINFO.longStart and (lon+dlon0)[0:npts/4-1] le flexINFO.longEnd)
      oplot,x[idx],y[idx] ,color=lineColor

 ;-Cross Horizontal Axis - Left
  idx=where((lon+dlon0)[npts/4:npts/2-1] ge flexINFO.longStart and (lon+dlon0)[npts/4:npts/2-1] le flexINFO.longEnd)+npts/4
      oplot,x[idx],y[idx] ,color=lineColor

 ;-Cross Vertical Axis 
  idx=where((lon+dlon0)[npts/2:*] ge flexINFO.longStart and (lon+dlon0)[npts/2-1:*] le flexINFO.longEnd)+npts/2
      oplot,x[idx],y[idx] ,color=lineColor






;--------------------;
;- Ellipse Plotting -;
;--------------------;
  for kPT=0,npts-1 do begin
     ux_arr[kPT]=1./cos(lat*!DTOR)*cos(2.*!DPI*kPT/float(npts))
     uy_arr[kPT]=1.*sin(2*!DPI*kPT/float(npts))
  endfor
  dtau=flexINFO.Radius/float(npts)
  dlat0=dlat0*0
  dlon0=dlon0*0
  for i=0,npts-1 do begin
     dlon0=dlon0+ux_arr*dtau
     dlat0=dlat0+uy_arr*dtau
     ux_arr=ux_arr+2*tan((lat+dlat0)*!DTOR)*!DTOR*ux_arr*uy_arr*dtau
     uy_arr=uy_arr-0.5*sin(2*(lat+dlat0)*!DTOR)*!DTOR*ux_arr*ux_arr*dtau
  endfor
          
  deproject,x,y,lat+dlat0,lon+dlon0,projection=projection,forward=1
  x = x mod (2.*!DPI)
  
  
  idx=where((lon+dlon0)ge flexINFO.longStart and (lon+dlon0)le flexINFO.longEnd)

  tmp=idx-shift(idx,1)
  idx2=where(tmp ne 1 and tmp ge 0,nbreak)
  if nbreak eq 0 then begin
     oplot,x[idx],y[idx],color=lineColor
  endif else begin
     oplot,x[idx[0:idx2[0]-1]],y[idx[0:idx2[0]-1]],color=lineColor
     oplot,x[idx[idx2[nbreak-1]:*]],y[idx[idx2[nbreak-1]:*]],color=lineColor
  endelse


endfor

end