clear all

%% Save locations Script
% Script can be rutned into a function later for simply ensuring filepath
% is valid for specific computer running matlab.
  MacFilePath='/Users/richarddudley/Dropbox/ParticleProductionDraft/';
    PCFilePath='D:\Dropbox\ParticleProductionDraft\';
    ClusterFilePath='/deac/generalGrp/andersonGrp/dudlra13/DensityDensity/ParticleProductionDraft/';
    
    
% If exist returns 7 if the path is a folder
    Mac=exist(MacFilePath,'dir')
    PC=exist(PCFilePath,'dir')
    Cluster=exist(ClusterFilePath,'dir')

    if Mac==7
        filepath = MacFilePath;
        filepath2= MacFilePath;

    elseif PC==7
        filepath = PCFilePath;

    elseif Cluster==7;
        filepath = ClusterFilePath;

    else
        filepath = pwd; % pwd is current working directory
    end


    
    
    %%This is for fixed asymtotic values of c1 and c2
    
%% ListMasses=[10^(-10), 10^(-6), 10^(-4), 10^(-3), 5*10^(-3), 10^(-2), 2*10^(-2) , 4*10^(-2), 6*10^(-2), 8*10^(-2), 10^(-1), 2*10^(-1), 3*10^(-1), 4*10^(-1)];% Mass list for data sets from PRD 2018

syms x omega VL VR
c1 = ((1/4)+(3/8).*(1+2.*pi.^(-1).*atan((1/8).*(8.*((-2)+3.^(1/2))+x)))).^(1/2);
kappa = double(subs(diff(c1, x), x, 0));


nuLZeroPotential=((-1)+exp(1).^(0.338502723697478160E3.*omega)).^(-1).*omega.^(-2) ...
  .*abs(omega).^2;

nuR2DeltaGeneric=4.*((-1)+exp(1).^(0.338502723697478160E3.*omega)).^(-1).*abs( ...
  omega.*((sqrt(-1)*2).*omega+(-1).*VR).^(-1)).^2;

nuL2DeltaGeneric=(1/4).*((-1)+exp(1).^(0.338502723697478160E3.*omega)).^(-1).* ...
  omega.^(-2).*(4.*omega.^2+VR.^2).^(-1).*(16.*omega.^4+(-1).*((-1)+ ...
  2.*exp(1).^(0.169251361848739080E3.*omega)+(-1).*exp(1).^( ...
  0.338502723697478160E3.*omega)).*VL.^2.*VR.^2+4.*omega.^2.*(exp(1) ...
  .^(0.338502723697478160E3.*omega).*VL.^2+2.*exp(1).^( ...
  0.169251361848739080E3.*omega).*VL.*VR+VR.^2));

    
nvL2DeltaGeneric=(1/4).*((-1)+exp(1).^(0.338502723697478160E3.*omega)).^(-1).* ...
  omega.^(-2).*(4.*omega.^2+VR.^2).^(-1).*((4.*omega.^2+VL.^2).* ...
  VR.^2+(-2).*exp(1).^(0.169251361848739080E3.*omega).*VL.*VR.*((-4) ...
  .*omega.^2+VL.*VR)+exp(1).^(0.338502723697478160E3.*omega).* ...
  VL.^2.*(4.*omega.^2+VR.^2));
    
   npoints=1:1/100:6;
   npoints2=1:1/100:7;


    
   %%
      HighOmega=10.^(npoints2/3-4/3);
    nuLZeroPotentialsubsHigha=subs(omega*nuLZeroPotential, omega,HighOmega);
    nuLZeroPotentialsubsHighb=subs(omega*nuLZeroPotential, omega,HighOmega)./nuLZeroPotentialsubsHigha(1);
       
         nuL2DeltaHigh1a=subs(nuL2DeltaGeneric,[VR,VL],[10*kappa, -100*kappa]);
     nuL2DeltaHigh1amax=subs(omega*nuL2DeltaHigh1a,omega,HighOmega);
    nuL2DeltaHigh1b=subs(omega*nuL2DeltaHigh1a,omega,HighOmega )./nuL2DeltaHigh1amax(1);
    
  %%  
 semilogx(HighOmega,nuLZeroPotentialsubsHighb*nuLZeroPotentialsubsHigha(1),'LineWidth',1,'Color', [0.6350 0.0780 0.1840])
 
 hold on   
  
semilogx( HighOmega,nuL2DeltaHigh1b*nuL2DeltaHigh1amax(1),'.-','LineWidth',1,'Color', [0 0.4470 0.7410])
 
 
 
hXLabel=xlabel('\omega'); 
    hYLabel=ylabel('C  \omega n_u^l '); 
    set([ hXLabel, hYLabel],'FontName', 'AvantGarde');
    set([ hXLabel, hYLabel], 'FontSize' , 16);
    ylim([0, 1.1])
    xlim([0.09, 11])
    graph1=gca;
    graph1.YTick=[0 0.5 1];
    graph1.YTickLabel = [{'0' ,'1/2', '1'}];
    
hold off
 print(gcf, '-djpeg', '-r500', [filepath,'/Figures/figure12a.png']);   
    
    
 
 
 semilogx(HighOmega,nuLZeroPotentialsubsHighb,'LineWidth',1,'Color', [0.6350 0.0780 0.1840])
 
 hold on   
  
semilogx( HighOmega,nuL2DeltaHigh1b,'.-','LineWidth',1,'Color', [0 0.4470 0.7410])
 
 
 
hXLabel=xlabel('\omega'); 
    hYLabel=ylabel('C  \omega n_u^l '); 
    set([ hXLabel, hYLabel],'FontName', 'AvantGarde');
    set([ hXLabel, hYLabel], 'FontSize' , 16);
    ylim([0, 1.1])
    xlim([0.09, 11])
    graph1=gca;
    graph1.YTick=[0 0.5 1];
    graph1.YTickLabel = [{'0' ,'1/2', '1'}];
    
hold off
 print(gcf, '-djpeg', '-r500', [filepath,'/Figures/figure12.png']);  
 
  loglog(HighOmega,nuLZeroPotentialsubsHighb,'LineWidth',1,'Color', [0.6350 0.0780 0.1840])
 
 hold on   
  
loglog( HighOmega,nuL2DeltaHigh1b,'.','LineWidth',1,'Color', [0 0.4470 0.7410])
 
 
 
hXLabel=xlabel('\omega'); 
    hYLabel=ylabel('C  \omega n_u^l '); 
    set([ hXLabel, hYLabel],'FontName', 'AvantGarde');
    set([ hXLabel, hYLabel], 'FontSize' , 16);
ylim([0.000001, 2])
%     xlim([0.09, 11])
%     graph1=gca;
%     graph1.YTick=[0 0.5 1];
%     graph1.YTickLabel = [{'0' ,'1/2', '1'}];
    
hold off
 print(gcf, '-djpeg', '-r500', [filepath,'/Figures/figure12V2.png']);  
 
 
 
   loglog(HighOmega,nuLZeroPotentialsubsHighb,'LineWidth',1,'Color', [0.6350 0.0780 0.1840])
 
 hold on   
  
loglog( HighOmega,nuL2DeltaHigh1b,'.','LineWidth',1,'Color', [0 0.4470 0.7410])
 
 
 
hXLabel=xlabel('\omega'); 
    hYLabel=ylabel('C  \omega n_u^l '); 
    set([ hXLabel, hYLabel],'FontName', 'AvantGarde');
    set([ hXLabel, hYLabel], 'FontSize' , 16);
ylim([10^-50, 20])
%     xlim([0.09, 11])
%     graph1=gca;
%     graph1.YTick=[0 0.5 1];
%     graph1.YTickLabel = [{'0' ,'1/2', '1'}];
    
hold off
 print(gcf, '-djpeg', '-r500', [filepath,'/Figures/figure12V3.png']);  
 
    loglog(HighOmega,nuLZeroPotentialsubsHighb*nuLZeroPotentialsubsHigha(1),'LineWidth',1,'Color', [0.6350 0.0780 0.1840])
 
 hold on   
  
loglog( HighOmega,nuL2DeltaHigh1b*nuL2DeltaHigh1amax(1),'.','LineWidth',1,'Color', [0 0.4470 0.7410])
 
 
 
hXLabel=xlabel('\omega'); 
    hYLabel=ylabel('\omega n_u^l '); 
    set([ hXLabel, hYLabel],'FontName', 'AvantGarde');
    set([ hXLabel, hYLabel], 'FontSize' , 16);
ylim([10^-50, 20])
%     xlim([0.09, 11])
%     graph1=gca;
%     graph1.YTick=[0 0.5 1];
%     graph1.YTickLabel = [{'0' ,'1/2', '1'}];
    
hold off
 print(gcf, '-djpeg', '-r500', [filepath,'/Figures/figure12V4.png']);  
 
 
 
    loglog(HighOmega,nuLZeroPotentialsubsHighb*nuLZeroPotentialsubsHigha(1),'LineWidth',1,'Color', [0.6350 0.0780 0.1840])
 
 hold on   
  
loglog( HighOmega,nuL2DeltaHigh1b*nuL2DeltaHigh1amax(1),'.','LineWidth',1,'Color', [0 0.4470 0.7410])
 
 
 
hXLabel=xlabel('\omega'); 
    hYLabel=ylabel('\omega n_u^l '); 
    set([ hXLabel, hYLabel],'FontName', 'AvantGarde');
    set([ hXLabel, hYLabel], 'FontSize' , 16);
ylim([10^-50, 20])
%     xlim([0.09, 11])
%     graph1=gca;
%     graph1.YTick=[0 0.5 1];
%     graph1.YTickLabel = [{'0' ,'1/2', '1'}];
    
hold off
 print(gcf, '-djpeg', '-r500', [filepath,'/Figures/figure12V4.png']);  
 
 %%   
    LowOmega=10.^(npoints-6);
     nuR2Delta1a=subs(nuR2DeltaGeneric,VR,10*kappa );
     nuR2Delta1amax=max(subs(omega*nuR2Delta1a,omega,LowOmega ));
    nuR2Delta1b=subs(omega*nuR2Delta1a,omega,LowOmega )./nuR2Delta1amax;
    
      nuR2Delta2a=subs(nuR2DeltaGeneric,VR,kappa );
      nuR2Delta2amax=max(subs(omega*nuR2Delta2a,omega,LowOmega ));
    nuR2Delta2b=subs(omega*nuR2Delta2a,omega,LowOmega )./nuR2Delta2amax;
    
      nuR2Delta3a=subs(nuR2DeltaGeneric,VR,2/3*kappa );
      nuR2Delta3amax=max(subs(omega*nuR2Delta3a,omega,LowOmega ));
    nuR2Delta3b=subs(omega*nuR2Delta3a,omega,LowOmega )./nuR2Delta3amax;
    
    %%
    
         nuL2Delta1a=subs(nuL2DeltaGeneric,[VR,VL],[10*kappa, -100*kappa] );
     nuL2Delta1amax=subs(omega*nuL2Delta1a,omega,LowOmega );
    nuL2Delta1b=subs(omega*nuL2Delta1a,omega,LowOmega )./nuL2Delta1amax(1);
    
      nuL2Delta2a=subs(nuL2DeltaGeneric,[VR,VL],[kappa, -100*kappa]);
      nuL2Delta2amax=subs(omega*nuL2Delta2a,omega,LowOmega );
    nuL2Delta2b=subs(omega*nuL2Delta2a,omega,LowOmega )./nuL2Delta2amax(1);
    
      nuL2Delta3a=subs(nuL2DeltaGeneric,[VR,VL],[2/3*kappa, -100*kappa]);
      nuL2Delta3amax=subs(omega*nuL2Delta3a,omega,LowOmega );
    nuL2Delta3b=subs(omega*nuL2Delta3a,omega,LowOmega )./nuL2Delta3amax(1);
    
    
    
  %%  
 semilogx(LowOmega,nuR2Delta3b,'LineWidth',1,'Color', [0.8500 0.3250 0.0980])
 
 hold on   
  
semilogx( LowOmega,nuR2Delta2b,'.','LineWidth',1,'Color', [0.4660 0.6740 0.1880])
 
semilogx( LowOmega,nuR2Delta1b,'--','LineWidth',1,'Color', [0 0.4470 0.7410])
 
 
hXLabel=xlabel('\omega'); 
    hYLabel=ylabel('C  \omega n_u^r '); 
    set([ hXLabel, hYLabel],'FontName', 'AvantGarde');
    set([ hXLabel, hYLabel], 'FontSize' , 16);
    ylim([0, 1.1])
    graph1=gca;
    graph1.YTick=[0 0.5 1];
    graph1.YTickLabel = [{'0' ,'1/2', '1'}];
    
hold off
 print(gcf, '-djpeg', '-r500', [filepath,'/Figures/figure13a.png']);
 
 

 %%   
 semilogx(LowOmega,nuL2Delta3b,'LineWidth',1,'Color', [0.8500 0.3250 0.0980])
 
 hold on   
  
semilogx( LowOmega,nuL2Delta2b,'.','LineWidth',1,'Color', [0.4660 0.6740 0.1880])
 
semilogx( LowOmega,nuL2Delta1b,'--','LineWidth',1,'Color', [0 0.4470 0.7410])
 
 
hXLabel=xlabel('\omega'); 
    hYLabel=ylabel('C  \omega n_u^l '); 
    set([ hXLabel, hYLabel],'FontName', 'AvantGarde');
    set([ hXLabel, hYLabel], 'FontSize' , 16);
    ylim([0, 1.1])
    graph1=gca;
    graph1.YTick=[0 0.5 1];
    graph1.YTickLabel = [{'0' ,'1/2', '1'}];
    
hold off
 print(gcf, '-djpeg', '-r500', [filepath,'/Figures/figure13b.png']);
    
    
    
    %%
    
     LowOmega2=0.008:0.0001:0.11;

     
     nuRmassive=4.*((-1)+exp(1).^(0.338502723697478160E3.*omega)).^(-1).*omega.*(( ...
  -7/10000)+omega.^2).^(1/2).*(omega+((-7/10000)+omega.^2).^(1/2)) ...
  .^(-2).*heaviside((-1/100).*7.^(1/2)+omega);
    
    nuLmassive=(1/4).*((-1)+exp(1).^(0.338502723697478160E3.*omega)).^(-1).* ...
  omega.^(-1).*((1/2000)+omega.^2).^(-1/2).*abs((-1).*omega+(( ...
  1/2000)+omega.^2).^(1/2)+exp(1).^(0.169251361848739080E3.*omega).* ...
  (omega+(-1).*((-7/10000)+omega.^2).^(1/2)).*(omega+((-7/10000)+ ...
  omega.^2).^(1/2)).^(-1).*((-1).*omega+((1/2000)+omega.^2).^(1/2))) ...
  .^2+((-7/10000)+omega.^2).^(1/2).*((1/2000)+omega.^2).^(-1/2).*( ...
  omega+((-7/10000)+omega.^2).^(1/2)).^(-2).*((-1).*omega+((1/2000)+ ...
  omega.^2).^(1/2)).^2.*heaviside((-1/100).*7.^(1/2)+omega);
    
    nuRmassivesubs=subs(nuRmassive, omega,LowOmega2 );
    
    nuLmassivesubs=subs(nuLmassive, omega,LowOmega );
    
    
    
    
    
    
    %%
    
     semilogx(LowOmega2,nuRmassivesubs,'LineWidth',1,'Color', [0.8500 0.3250 0.0980])
 
 
hXLabel=xlabel('\omega'); 
    hYLabel=ylabel(' n_u^r '); 
    set([ hXLabel, hYLabel],'FontName', 'AvantGarde');
    set([ hXLabel, hYLabel], 'FontSize' , 16);
xlim([0.009, 0.11])
%     graph1=gca;
%     graph1.YTick=[0 0.5 1];
%     graph1.YTickLabel = [{'0' ,'1/2', '1'}];
    
 print(gcf, '-djpeg', '-r500', [filepath,'/Figures/figure18a.png']);
    
    
    
    
    
    
         semilogx(LowOmega,nuLmassivesubs,'LineWidth',1,'Color', [0.8500 0.3250 0.0980])
 
 
hXLabel=xlabel('\omega'); 
    hYLabel=ylabel('n_u^l '); 
    set([ hXLabel, hYLabel],'FontName', 'AvantGarde');
    set([ hXLabel, hYLabel], 'FontSize' , 16);
%     graph1=gca;
%     graph1.YTick=[0 0.5 1];
%     graph1.YTickLabel = [{'0' ,'1/2', '1'}];
    
 print(gcf, '-djpeg', '-r500', [filepath,'/Figures/figure18b.png']);
    
    
  %%
  VNaught=3/4;
  x3=10;
  options = odeset('RelTol',2.23e-12,'AbsTol', 1e-12); 
  points=200;
  xstarOut=zeros(points,1);
  xstarIn=zeros(points,1);
  xin=zeros(points,1);
  xout=zeros(points,1);

 
  for i=1:points
      xmax=i/10.00003;
                x1span=[x3,xmax]; 
                x1spanIn=[-x3,-xmax]; % span for ode value, last value in array is equal to the value desired.
                Y10=0;  % Initial value for tortise coordinate
                xStar1Solver = @(x,y) (sqrt((1/2)^2+1/2*((1)^2-(1/2)^2)*(1+(2/pi)*atan((x+(8*tan((pi/((1)^2-(1/2)^2))*((3/4)^2-1/2*((1/2)^2)+1^2))))/(8)))))./((1/2)^2+1/2*((1)^2-(1/2)^2)*(1+(2/pi)*atan((x+(8*tan((pi/((1)^2-(1/2)^2))*((3/4)^2-1/2*((1/2)^2)+1^2))))/(8)))-(VNaught)^2);
                [T1,Y1] = ode45(xStar1Solver,x1span,Y10,options);
                xstarOut(i)=Y1(end);
                
                xStar1SolverIn = @(x,y) (sqrt((1/2)^2+1/2*((1)^2-(1/2)^2)*(1+(2/pi)*atan((x+(8*tan((pi/((1)^2-(1/2)^2))*((3/4)^2-1/2*((1/2)^2)+1^2))))/(8)))))./((1/2)^2+1/2*((1)^2-(1/2)^2)*(1+(2/pi)*atan((x+(8*tan((pi/((1)^2-(1/2)^2))*((3/4)^2-1/2*((1/2)^2)+1^2))))/(8)))-(VNaught)^2);
                [T1in,Y1in] = ode45(xStar1SolverIn,x1spanIn,Y10,options);
                xstarIn(i)=Y1in(end);
        xin(i)=(-xmax);
        xout(i)=(xmax);
  end
    xtotal = cat(1, flip(xin), xout);
    xstartotal=cat(1, flip(xstarIn), xstarOut);
    LineSize=size(xstartotal);
    horizontalLine=zeros(LineSize(1), 1);
    
     plot(xtotal,xstartotal,'LineWidth',1,'Color', [0.8500 0.3250 0.0980])
     hold on
     plot(xtotal,horizontalLine,'LineWidth',1,'Color', [0.0 0.0 0.0])
line(zeros(11,1), -100:20:100);
 txt0 = 'Horizon';
 text(-4.5,35,txt0)
 txt1 = 'x_r';
text(9.75,-7,txt1)
 txt1a = '|';
text(10,0,txt1a)
txt2 = 'x_l';
text(-10.25,-7,txt2)
text(-10,0,txt1a)
 
hXLabel=xlabel('x'); 
    hYLabel=ylabel('x^*'); 
    set([ hXLabel, hYLabel],'FontName', 'AvantGarde');
    set([ hXLabel, hYLabel], 'FontSize' , 16);
ylim([-100, 40])
%     graph1=gca;
%     graph1.YTick=[0 0.5 1];
%     graph1.YTickLabel = [{'0' ,'1/2', '1'}];
    
hold off
 print(gcf, '-djpeg', '-r500', [filepath,'/Figures/figure3.png']);
    