%function plot_Earth(opt,Ls,MLTs)
function plot_Earth(plane,varargin)
% example here
% figure; plot_Earth('2M','l',[3,4,5]); xlim([0,6])
% Lunjin Chen
if nargin==0,plane='2E';,end
B0=0.31d-4; % Earth dipole
Ls=[];
Ld0s=[];
Ls_linedip=[];
dB=0;
MLTs=[];
Rs=[];
LATs=[];
i = 1;
latlnsty=':k';
while i <= length(varargin)
    switch lower(varargin{i})
        case 'l'
            Ls  =varargin{i+1};   i=i+2;
        case 'l_linedip'
            Ls_linedip  =varargin{i+1};   i=i+2;            
        case 'db'
            dB  =varargin{i+1};   i=i+2;
        case 'ld0'
            Ld0s  =varargin{i+1};   i=i+2;
        case 'mlt'
            MLTs  =varargin{i+1};   i=i+2;
        case 'r'
            Rs  =varargin{i+1};   i=i+2;
        case 'lat'
            LATs  =varargin{i+1};   i=i+2;
        case 'latlnsty'
            latlnsty=varargin{i+1};   i=i+2;
        otherwise
            error(['Input argument ',varargin{i},' not valid.']);
    end
end


if (strexist_(plane,'BT')) %magnetic tube
    x=[];
    z=[];
    nL=length(Ls);
    for iL=1:nL
        L=Ls(iL);
        ILAT=acos(sqrt(1./L));
        lat=linspace(-ILAT,ILAT,100);
        r=L.*cos(lat).^2;
        x{iL}=r.*cos(lat);
        z{iL}=r.*sin(lat);
    end
    if (nL==1)
        x{2}=x{1}([end,1]);
        z{2}=z{1}([end,1]);
    end
    clcarr=gray(3);
    patch([x{1},fliplr(x{2})],[z{1},fliplr(z{2})],clcarr(2,:),'edgecolor','none','facealpha',0.4)
    
end

if (strexist_(plane,'2M')) % 2D
    hold on
    ylabel('Z, R_E'); xlabel('R_{XY}, R_E')
    for L=Ls
        ILAT=acos(sqrt(1./L));
        lat=linspace(-ILAT,ILAT,100);
        r=L.*cos(lat).^2;
        x=r.*cos(lat);
        z=r.*sin(lat);
        plot(x,z,'--k')
        %plot(x,z,'-k','linewidth',2)
    end
    
    for L=Ls_linedip
        ILAT=acos(1./L);
        lat=linspace(-ILAT,ILAT,100);
        r=L.*cos(lat);
        x=r.*cos(lat);
        z=r.*sin(lat);
        plot(x,z,'--k')
        %plot(x,z,'-k','linewidth',2)
    end
    
    for L=Ld0s
        bfac=dB./B0;
        % what happens if dB=-B0./L^3
        C=(1-bfac/2.*L.^3)./L;
        if (bfac*L.^3>=-1)
            r=linspace(1,L,1000);
        else
            r=linspace(L,2*L,1000);
        end
        clat2=C.*r./(1-bfac/2.*r.^3);
        clat=clat2*NaN;
        ind=clat2>=0 & clat2<=1;
        clat(ind)=sqrt(clat2(ind));
        lat=acos(clat);
        x=r.*cos(lat);
        z=r.*sin(lat);
        plot([fliplr(x),x],[-fliplr(z),z],'--k')
        %plot(x,z,'-k','linewidth',2)
    end
    
    for lat=LATs
        x=[1.05 12].*cos(lat*pi/180);
        z=[1.05 12].*sin(lat*pi/180);
        plot(x,z,latlnsty)
        %plot(x,-z,':k')
    end
    
    if ~(strexist_(plane,'_noearth_'))    
    theta1=linspace(-pi/2,pi/2,100);
    theta2=linspace(pi/2,pi*3/2,100);
    
    fill(cos(theta1),sin(theta1),'k')
    plot(cos(theta2),sin(theta2),'-k')
    end
    axis equal
    hold off
end

if (strexist_(plane,'Earth')) % 2D Dark Earth
    hold on
    theta1=linspace(-pi/2,pi/2,100);
    theta2=linspace(pi/2,pi*3/2,100);
    fill(cos(theta2),sin(theta2),'-k')
    fill(cos(theta1),sin(theta1),'k')
    axis equal
    hold off
end


if (strexist_(plane,'2E')) % 2D
    hold on
    theta1=linspace(-pi/2,pi/2,100);
    theta2=linspace(pi/2,pi*3/2,100);
    plot(cos(theta2),sin(theta2),'-k')
    fill(cos(theta1),sin(theta1),'k')
    for MLT=MLTs
        x=[1.01 15]*cos(MLT*pi/12);
        y=[1.01 15]*sin(MLT*pi/12);
        plot(x,y,'--k')
    end
    for r=Rs
        theta3=linspace(-pi,pi,100);
        x=r*cos(theta3);
        y=r*sin(theta3);
        plot(x,y,'--k')
    end
    
    axis equal
    hold off
end


if (strexist_(plane,'3D')) % 3D
    hold on
    sphere(30)
    for L=Ls
        for MLT=MLTs
            ILAT=acos(sqrt(1./L));
            lat=linspace(-ILAT,ILAT,100);
            r=L.*cos(lat).^2;
            x=r.*cos(lat).*cos(MLT*pi/12);
            y=r.*cos(lat).*sin(MLT*pi/12);
            z=r.*sin(lat);
            plot3(x,y,z,'--k')
        end
    end
    axis equal
    hold off
end


if (strexist_(plane,'fill'))
    theta1=linspace(-pi,pi,100);
    fill(cos(theta1),sin(theta1),'k')
end

if (strexist_(plane,'lfill'))
    for L=Ls(1)
        ILAT=acos(sqrt(1./L));
        lat=linspace(-ILAT,ILAT,100);
        r=L.*cos(lat).^2;
        x=r.*cos(lat);
        z=r.*sin(lat);
        %plot(x,z,'--k')
        %plot(x,z,'-k','linewidth',2)
        %clcarr=colormap(gray(3))
        clcarr=gray(3)
        fill(x,z,clcarr(2,:))
    end
    
end



function flag=strexist_(str,substr)
% check whether substr exists in str
% 12/20/2013 Lunjin Chen
flag=~isempty(strfind(str,substr));
return

