function ostr=tpgcdip(pstr,wstr,flag)
% test particle averaged over gyro motion	  
% field-aligned waves, dipole field 
% input/output angle in angle
% angle during numerical calculation in unit of radian
% Lunjin Chen 1/14/2013
global latlast_
ostr=[];

if nargin==0
    pstr.iq=-1;
    pstr.imas=0;
    pstr.L0=5;
    pstr.lat0=0;
    pstr.latlast=25;
    pstr.en0=1e6; % kinetic energy in eV
    pstr.pa0=4; % pitch angle in deg
    pstr.ga0=0; %  gyro angle in deg
    
    pstr.t=linspace(0,0.01,10000); % sec

    wstr.ieta=[0.8 0.2];
    wstr.ne_func=@(lat) 50d6./cos(lat).^4; % Denton et al., 2002
    wstr.ne_func=@(lat) 50d6;
    wstr.freq=3.9/16;
    wstr.lrmode='L';
    wstr.bw_func=@(lat) 2d-9 ; % wave magnetic field
end


phyval=calphy();

abstol=1d-9;
reltol=1d-9;
ev=phyval.ev;
c=phyval.c;

q=pstr.iq*ev;
if pstr.imas==0
	mas=phyval.emas;
else
	mas=phyval.hmas*pstr.imas;
end

if nargin<3
	flag=1;
end

if flag==0
	flagdip=1;
	flaguni=0;
else	
	flagdip=0;
	flaguni=1;
end

L0=pstr.L0; % initial L
lat0=pstr.lat0*pi/180; % initial laitude
latlast_=pstr.latlast*pi/180; % final latitude

gam0=pstr.en0*ev/(mas*c^2)+1;
ga0=pstr.ga0*pi/180;
pa0=pstr.pa0*pi/180;
Vt0=c*sqrt(1-1/gam0^2);
Vx0=Vt0*sin(pa0)*cos(ga0);
Vy0=Vt0*sin(pa0)*sin(ga0);
Vz0=Vt0*cos(pa0);
Vp0=Vt0*sin(pa0);

Px0=gam0*mas*Vx0;
Py0=gam0*mas*Vy0;
Pz0=gam0*mas*Vz0;
Pp0=gam0*mas*Vp0;
Pt0=gam0*mas*Vt0;

% z: //B0, p: perpendicular to B0, t: total
% eq: equatorial value


B0eq=phyval.B0dip/L0^3;
R0=phyval.RE*L0;
Wc0=abs(q)*B0eq/mas;

mc=mas*phyval.c;
mc2=mc^2;
if strcmpi(wstr.lrmode,'L')
	sgnlr=1; % L-mode
else
	sgnlr=-1; % for R-mode
end



t=pstr.t;

Rn=R0;
Pn=Pt0;

etaH=wstr.ieta(1);
etaHe=wstr.ieta(2);
etaO=1-etaH-etaHe;
ne_func_=wstr.ne_func; % plasma density function
bw_func_=wstr.bw_func; % wave magnetic field 
b0_func_=@(lat)B0eq*sqrt(1+3*sin(lat)^2)/cos(lat)^6; % dipole field
b0_func_=@(lat)B0eq;


phyval=calphy(struct('ne',ne_func_(0),'b',B0eq));
if(wstr.freq>0)
	  freq=phyval.fch*wstr.freq; % normalized frequency
else
	  freq=abs(wstr.freq); % absolute frequency
end
w=2*pi*freq;
%kz_func_=@(lat)kfaw_(w,ntot1,B0,etaH,etaHe,etaO,sgnlr);
[w,ne_func_(0),b0_func_(0),etaH,etaHe,etaO,sgnlr]
kz_func_=@(lat)kfaw(w,ne_func_(lat),b0_func_(lat),etaH,etaHe,etaO,sgnlr);


% setup initial condition
eta0=ga0;
z0=z_dip_(lat0);
y0=[Pz0,Pp0,eta0,z0,lat0];

fevents=@events_latlast;
opt=odeset('AbsTol',abstol*[Pn,Pn,1,Rn,1],'RelTol',reltol,'Events',fevents);
[to,yo]=ode45(@tpode_,t,y0,opt);

ostr=pstr;
ostr.wstr=wstr;
ostr.ne_func=ne_func_;
ostr.bw_func=bw_func_;
ostr.b0_func=b0_func_;
ostr.kz_func=kz_func_;
ostr.w=w;
ostr.Wc0=Wc0;


ostr.Vt0=Vt0;
ostr.Vp0=Vp0;
ostr.Vz0=Vz0;
ostr.gam0=gam0;
ostr.Pp0=Pp0;
ostr.Pz0=Pz0;
ostr.Pt0=Pt0;
ostr.t=to;
ostr.pz=yo(:,1);
ostr.pp=yo(:,2);
pt2=ostr.pz.^2+ostr.pp.^2;
ostr.pt=sqrt(pt2);
ostr.gam=sqrt(pt2./mc2+1);
ostr.pa=atan2(ostr.pp,ostr.pz)*180/pi;
ostr.eta=yo(:,3)*180/pi;
ostr.z=yo(:,4);
ostr.lat=yo(:,5)*180/pi;

plot_ostr_(ostr)

function yp=tpode_(t,y)
	  % [pz, pp, eta, z, lat]
	  % 
	  pz1=y(1);
	  pp1=y(2);
	  eta1=y(3);
	  z1=y(4);
	  lat1=y(5);

	  pp2=pp1.^2;
	  pz2=pz1.^2;
	  p2=pz2+pp2;
	  g1=sqrt(p2./mc2+1);
	  gmas1=g1*mas;
	  vz1=pz1/gmas1;

	  seta1=sin(eta1);
	  ceta1=cos(eta1);

	  slat1=sin(lat1);
	  clat1=cos(lat1);

	  termlat=sqrt(1+3*slat1^2);


	  % dzdlat=dzfunc_(lat1);
	  dzdlat=R0*clat1*termlat;

	  if flagdip
	  	%B0=b0func_(lat1);
	  	B0=B0eq*termlat/clat1^6; % dipole field
	  	%db0term=db0func_(lat1);
	  	db0dlat=B0eq*3*slat1*(3+5*slat1^2)/(clat1^7*termlat); % dipole field
	  	db0term=db0dlat/dzdlat/B0; % dipole field
	  elseif flaguni % uniform case
	  	B0=B0eq;
		db0term=0.0d0;
	  else
	  end

	  %kz1=kz_func_(lat1);
	  ntot1=ne_func_(lat1); % plasma density
	  %kz1=kfaw_(w,ntot1,B0,etaH,etaHe,etaO,sgnlr);
	  % kz1>0, wave moving in positive z direction
	  kz1=kz_func_(lat1);


	  % field-aligned waves
	  Bw1=bw_func_(lat1);
	  Ew1=w/kz1*Bw1;

	  qBw1=q*Bw1/gmas1;
	  qEw1=q*Ew1;
	  yp=zeros(5,1);
	  yp(1)=-qBw1*pp1*seta1-(pp2/gmas1*0.5)*db0term;
	  yp(2)=(-qEw1+qBw1*pz1)*seta1+(pz1*pp1/gmas1*0.5)*db0term;
	  yp(3)=(-qEw1/pp1+qBw1*pz1/pp1)*ceta1-q*B0/gmas1+  sgnlr*(kz1*vz1-w);
	  yp(4)=vz1;
	  yp(5)=vz1/dzdlat;
end	  


end	  

function B0=B0fun_dip_(lat)
	  % variation of B0 normalized to B0eq
end	  

function z=z_dip_(lat)
	  % z normalized to R0
	  t1=sin(lat);
	  z=0.5*t1.*sqrt(3*t1.^2+1)+sqrt(3)/6*asinh(sqrt(3)*t1);
end	  

function k=kfaw_(w,ntot,btot,eta1,eta2,eta3,sgnlr)
hmas=1.67252e-27;
ev  =1.60210e-19;
emas=9.10908e-31;
e0  =8.85442e-12;
ee0 =e0;
mu0 =4.0*pi*1e-7;
clight=2.99792501e8;
kb  =1.38054e-23;
rad =pi/180.;
epsilon=emas/hmas;
memp=epsilon;


wpe=sqrt(ev*ev/emas/e0*ntot);
wce=ev/emas*btot;
wch=ev/hmas*btot;
fpara=wpe./wce;
x=w./wch;

alphastar=(1./fpara).^2; % 1 X nx

if sgnlr==1
	L=1.0-1./(alphastar.*memp.*x).* ...
    (1./(memp.*x+1)+eta1./(x-1)+eta2./(4*x-1)+eta3./(16*x-1)); % check
	k=sqrt(L).*w/clight;    
else    
	R=1.0-1./(alphastar.*memp.*x).* ...
    (1./(memp.*x-1)+eta1./(x+1)+eta2./(4*x+1)+eta3./(16*x+1));
	k=sqrt(R).*w/clight;    
end

end    

function [value,isterminal,direction] = events_latlast(t,y)
	% [value,isterminal,direction] = events(t,y). For the ith event function in events,
	% value(i) is the value of the function.
	% isterminal(i) = 1, if the integration is to terminate at a zero of this event function i
	% and 0 otherwise.
	% direction(i) = 0 if all zeros are to be computed (the default), 
	%   	          +1 if only the zeros where the event function increases, 
	%        	  and -1 if only the zeros where the event function decreases.
	global latlast_
    % slatlast_ is the maximum latitudes
	value = [y(5)-latlast_];     % detect r=1
	isterminal = [1];   % stop the integration
	direction = [0];   % all direction
end


function plot_ostr_(ostr)
	  figure(1)
	  clf
	  setfigpos(gcf,[1 1 4 10])
	  subplot(511)
	  plot(ostr.t,ostr.gam)
	  subplot(512)
	  plot(ostr.t,ostr.pa)
	  subplot(513)
	  plot(ostr.t,ostr.lat)
	  subplot(514)
	  plot(ostr.t,ostr.z)
	  subplot(515)
	  plot(ostr.t,ostr.eta)
	  saveeps('_')
end

