;+
; NAME:
;     LMN_TRANSFORM_CRIB
;
; PURPOSE:
;	This code shows anatomy of the LMN transform, so essentially this is dissected gsm2lmn
;	routine with hand-picked solar wind input. It projects the LMN base vectors onto the XY (GSM)
;	plane at the points of magnetopause boundary for different levels of Z (smooth line - Z=0,
;	dotted - Z=-10, dashed - Z=15 (Earth radii). The LMN orts are plotted in the sequence N-M-L,
;	so, if you see that the green one overshadows the red, that means that the green vector 
;	screens part of the red vector (you look from atop of the Z axis).
;
; CATEGORY:
;	Crib
;
; CALLING SEQUENCE:
;	lmn_transform_crib
;
; INPUTS:
;	none
;
; KEYWORDS: none
;
; PARAMETERS: none; image parameters are set in the code
;
; OUTPUTS:
;	graphic file mpfig.ps (color postscript)
;
; DEPENDENCIES: CDAWlib files DeviceOpen.pro and DeviceClose.pro.
;
; MODIFICATION HISTORY:
;	Written by:	Vladimir Kondratovich 2008/01/30
;-
;
; THE CODE BEGINS:


function equ,t
  compile_opt hidden
  common eqrhs,rhs,alpha
  const=exp(alpha*alog(2.))
  diff=const*sin(t)/exp(alpha*alog(1.+cos(t)))-rhs
  return,diff
end



PRO lmn_transform_crib

common eqrhs,rhs,alpha

;SW input:
dp=2.088; nPa
bz=0.; nT

;Shue et al. (1998) MP:
alpha = (0.58-0.007*bz)*(1+0.024*alog(dp))
r0=(10.22+1.29*tanh(0.184*(bz+8.14)))/exp(alog(dp)/6.6)

;for fun
yarr1=1.-2.*findgen(11)/10.
xarr1=sqrt(1.-yarr1^2)
yarr2=reverse(yarr1)
xarr2=fltarr(11)
yarr=[yarr1,yarr2]
xarr=[xarr1,xarr2]

;number of points to plot
nwing=14
npts=2*nwing+1

;we will plot 3 z-levels:
zarr=[0,-10,15]
linearr=[0,1,2]
nlines=n_elements(zarr)

;look for this file in your work directory after lmn_transform_crib finishes
nam='mpfig.ps'
deviceopen,1,/portrait,fileoutput=nam

red = [0,1,1,0,0,1]
green = [0,1,0,1,0,1]
blue = [0,1,0,0,1,0]
tvlct, 255 * red, 255 * green, 255 * blue

for j=0,nlines-1 do begin;+++++++++++++++++++++++++++++++++++++++++++++++++++++

  ;calculate X, Y, Z in GSM
  z0j=zarr(j)
  rhs=abs(z0j/r0)
  if rhs le 0.9 then t0=asin(rhs) else t0=1.5
  tt=newton(t0,'equ')
  th0=tt
  
  th=(findgen(nwing)+1)*!pi/18.+th0
  r1=r0*exp(alpha*alog(2./(1+cos(th))))
  ro=r0*exp(alpha*alog(2./(1+cos(th0))))
  r=[reverse(r1),ro,r1]
  
  x=fltarr(npts)
  x(nwing)=ro*cos(th0)
  x(0:nwing-1)=reverse(r1*cos(th))
  x(nwing+1:npts-1)=r1*cos(th)
  
  z=fltarr(npts)+z0j
  
  argu=r^2-z^2-x^2>0.
  y=sqrt(argu)
  y(0:nwing-1)=-y(0:nwing-1)
  
  vl=fltarr(npts,3)
  vm=fltarr(npts,3)
  vn=fltarr(npts,3)
  
  ;Transformation starts
  for i = 0, npts-1 do begin
     xi=x(i)
     yi=y(i)
     zi=z(i)
  
     theta=acos(xi/sqrt(xi^2+yi^2+zi^2));
     rho=sqrt(yi^2+zi^2)
     if rho gt 0. then begin
        tang1 = [0.,zi,-yi]
        tang2 = [xi,yi,zi]*alpha*sin(theta)/(1+cos(theta))+[-rho^2,xi*yi,xi*zi]/rho
        tang1 = tang1/sqrt(total(tang1^2));
        tang2 = tang2/sqrt(total(tang2^2));
        dN = crossp(tang1,tang2);
        dM = crossp(dN, [0,0,1])
        dM = dM/sqrt(total(dM^2))
        dL = crossp(dM, dN)
     endif else begin
        dN=[1.,0.,0.]
        dM=[0.,-1.,0.]
        dL=[0.,0.,1.]
     endelse
     transm = [[dL],[dM],[dN]]
  
     vl(i,*)=dL
     vm(i,*)=dM
     vn(i,*)=dN
  endfor
  
  ;the length of displayed LMN orts
  mult=1.
  xmin=min(x,max=xmax)
  ymin=min(y,max=ymax)
  
  if j eq 0 then begin
     plot, x, y, thick=1, charsize=1.3,line=linearr(j),$
        xrange=[xmin-mult, xmax+mult], xstyle = 8, ystyle = 8, $
        yrange=[ymin-mult, ymax+mult],$
        xtitle='X, GSM',ytitle='Y,  GSM',$
        title='LMN triades, L=green, M=blue, N=red',/noerase
     axis,0,0,xax=0,/data
     axis,0,0,0,yax=0,/data
     plots,xarr,yarr
     polyfill,xarr,yarr,color=1
     plots,-xarr,yarr
     polyfill,-xarr,yarr,color=0
  endif else begin
     oplot, x, y, thick=1,line=linearr(j)
  endelse
  
  for i=0,npts-1 do begin
     oplot,[x(i),x(i)+mult*vn(i,0)],[y(i),y(i)+mult*vn(i,1)],linestyle=0,color=2,thick = 3
     oplot,[x(i),x(i)+mult*vm(i,0)],[y(i),y(i)+mult*vm(i,1)],linestyle=0,color=4,thick = 3
     oplot,[x(i),x(i)+mult*vl(i,0)],[y(i),y(i)+mult*vl(i,1)],linestyle=0,color=3,thick = 3
  endfor

endfor;------------------------------------------------------------------------

deviceclose

print, 'lmn_transform_crib finished.'

end