%function f = DDToyHLower(TypePH, Region1PH, Region2PH, MassPH)
% Region 41 is the test case
for snark=1:1
for qart=2:2
    
TypePH='Large';
Region1PH=4;
Region2PH=1;
MassPH=2*qart*10^(-(7-snark));
%%Medium run at MassPH=4*10^(-2) takes 1.5 hours
%%Medium run at MassPH=4*10^(-3) takes 1 hours
%% 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/BH/BEC/Massive/DensityDensity/AnalyticModeSolutions/ToyModel/Total/';
    PCFilePath='D:\Dropbox\BH\BEC\Massive\DensityDensity\AnalyticModeSolutions\ToyModel\Total\';
    ClusterFilePath='/deac/generalGrp/andersonGrp/dudlra13/DensityDensity/InputData/';
    
    
% 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= '/Users/richarddudley/Dropbox/BH/BEC/Massive/DensityDensity/InputData/';
        addpath('/Users/richarddudley/Dropbox/BH/BEC/Massive/DensityDensity/AnalyticModeSolutions/ToyModel/AnalyticChecks/DD Function/ModeDefinitions/');
        addpath('/Users/richarddudley/Dropbox/BH/BEC/Massive/');
        savefilepath='/Users/richarddudley/Dropbox/BH/BEC/Massive/DensityDensity/TempOutput/';
    elseif PC==7
        filepath = PCFilePath;
        filepath2 = 'D:\Dropbox\BH\BEC\Massive\DensityDensity\InputData\';
        addpath('D:\Dropbox\BH\BEC\Massive\');
        addpath('D:\Dropbox\BH\BEC\Massive\DensityDensity\AnalyticModeSolutions\ToyModel\AnalyticChecks\DD Function\ModeDefinitions\');
        savefilepath='D:\Dropbox\BH\BEC\Massive\DensityDensity\TempOutput\';
    elseif Cluster==7;
        filepath = ClusterFilePath;
        addpath('/deac/generalGrp/andersonGrp/dudlra13/Functions/') 
    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
mass=MassPH;%sets mass value
mR = mass*sqrt(7/16);  % Effective mass in exterior region
mL = mass*sqrt(5/16); 
massString=num2str(mass);

RegionPHPH=strjoin({num2str(Region1PH),num2str(Region2PH)});
RegionPHPHPH= RegionPHPH(~isspace(RegionPHPH));
RegionPH=str2num(RegionPHPHPH);
RegionString=num2str(RegionPH);
% The following defintions are set for xstar and t calculations for input
% and must be adjusted if one wishes to modify the speed of sound.
syms x;
vNaught=3/4;
c1=1;
c2=1/2;
sigmaV=8;
b=sigmaV*(tan((pi/(c2^2-c1^2))*(vNaught^2-1/2*(c1^2+c2^2))));
c=sqrt(c1^2+1/2*(c2^2-c1^2)*(1+2/pi*atan((x+b)/sigmaV)));
dcdx=(c2^2-c1^2)/(2*pi*sigmaV*(1+(x+b)^2/sigmaV^2)*c);
kappaPH=subs(dcdx, x,0);
kappa=eval(kappaPH);% Represents the surface 'gravity' at the horizon
Lambda=1;



%Importing based on Region and type
xstarPH=importdata([filepath2, 'XandT/X',num2str(Region1PH),TypePH,'.txt']);
xstarPrimePH=importdata([filepath2, 'XandT/X',num2str(Region2PH),'Prime',TypePH,'.txt']);

    % Settings for ODE solvers   
        options = odeset('RelTol',2.23e-12,'AbsTol', 1e-12); 





IRcutoff=mR*0.00000000010; % Infrared cutoff

[nn, Trash]=size(xstarPH(:,1));
[mm, Trash]=size(xstarPrimePH(:,1));

mmm=mm;
nnn=nn;

AAAAA=zeros(nnn,mmm);
AAAAALR=zeros(nnn,mmm);
AAAAA2=zeros(nnn,mmm);
AAAAALR2=zeros(nnn,mmm);
        
    for n=1:nnn
        for m=1:mmm
            x1=xstarPH(n,1);
            xstar1=x1;
            xActual1=xstarPH(n,3);
            x2=xstarPrimePH(m,1);
            xstar2=x2;
            xActual2=xstarPrimePH(n,3);
            t1=xstarPH(n,2);
            t2=xstarPrimePH(m,2);
            c1PH=subs(c, x, x1);
            c1=eval(c1PH);
            
        if RegionPH==11
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable11LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable11Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable11LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable11Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==12        
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable12LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable12Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable12LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable12Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==13        
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable13LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable13Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable13LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable13Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==14       
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable14LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable14Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable14LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable14Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==21        
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable21LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable21Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable21LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable21Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==22
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable22LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable22Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable22LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable22Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==23       
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable23LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable23Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable23LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable23Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==24        
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable24LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable24Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable24LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable24Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==31       
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable31LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable31Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable31LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable31Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==32       
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable32LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable32Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable32LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable32Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==33        
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable33LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable33Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable33LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable33Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==34        
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable34LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable34Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable34LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable34Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==41       
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable41LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable41Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable41LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable41Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==42        
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable42LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable42Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable42LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable42Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==43        
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable43LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable43Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable43LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable43Upper(x1, x2, t1, t2, mR, mL, kappa);

        elseif RegionPH==44        
 %          [ILRLower, ILRttLower, ILRxxLower, ILRxtLower, ILRtxLower, ILRx1Lower, ILRx2Lower, ILRt1Lower, ILRt2Lower]=IDataTable44LRLower (x1, x2, t1, t2, mR, mL, kappa);
 %          [ILower, IttLower, IxxLower, IxtLower, ItxLower, Ix1Lower, Ix2Lower, It1Lower, It2Lower]=IDataTable44Lower (x1, x2, t1, t2, mR, mL, kappa);
              [ILRUpper, ILRttUpper, ILRxxUpper, ILRxtUpper, ILRtxUpper, ILRx1Upper, ILRx2Upper, ILRt1Upper, ILRt2Upper]=IDataTable44LRUpper(x1, x2, t1, t2, mR, mL, kappa);
              [IUpper, IttUpper, IxxUpper, IxtUpper, ItxUpper, Ix1Upper, Ix2Upper, It1Upper, It2Upper]=IDataTable44Upper(x1, x2, t1, t2, mR, mL, kappa);
        else
            
        end
   
        
            tspan=[mR,Lambda];
            I0=0;


                ISolver=IUpper;              
                [Omega,ILower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(IttUpper);              
                [Omega,IttLower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(IxxUpper);              
                [Omega,IxxLower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(IxtUpper);              
                [Omega,IxtLower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(ItxUpper);              
                [Omega,ItxLower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(Ix1Upper);              
                [Omega,Ix1Lower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(Ix2Upper);              
                [Omega,Ix2Lower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(It1Upper);              
                [Omega,It1Lower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(It2Upper);              
                [Omega,It2Lower1] = ode45(ISolver,tspan,I0, options);
                
                
            
                ISolver=ILRUpper;              
                [Omega,ILRLower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(ILRttUpper);              
                [Omega,ILRttLower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(ILRxxUpper);              
                [Omega,ILRxxLower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(ILRxtUpper);              
                [Omega,ILRxtLower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(ILRtxUpper);              
                [Omega,ILRtxLower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(ILRx1Upper);              
                [Omega,ILRx1Lower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(ILRx2Upper);              
                [Omega,ILRx2Lower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(ILRt1Upper);              
                [Omega,ILRt1Lower1] = ode45(ISolver,tspan,I0, options);
                
                ISolver=(ILRt2Upper);              
                [Omega,ILRt2Lower1] = ode45(ISolver,tspan,I0, options);
                
% here the coeficeints (hbar n)/(2m l) are neglected
c=sqrt(c1^2+1/2*(c2^2-c1^2)*(1+2/pi*atan((x+b)/sigmaV)));
dcdx=(c2^2-c1^2)/(2*pi*sigmaV*(1+(x+b)^2/sigmaV^2)*c);

cofx1=subs(c, x,xActual1);
cofx2=subs(c, x,xActual2);
dcdx1=subs(dcdx, x,xActual1);
dcdx2=subs(dcdx, x,xActual2);
cofxsquaredterm=1./((cofx1.^2).*(cofx1.^2));
dxstar1dx=cofx1/(cofx1.^2-vNaught^2);
dxstar2dx=cofx2/(cofx2.^2-vNaught^2);

t1x2coeff=(vNaught)*dxstar2dx*sqrt(cofx1*cofx2);
x1t2coeff=(vNaught)*dxstar1dx*sqrt(cofx1*cofx2);
x1x2coeff=(vNaught^2)*dxstar2dx*dxstar1dx*sqrt(cofx1*cofx2);
t1coeff=(vNaught)*dcdx2*sqrt(cofx1)/(2*sqrt(cofx2));
t2coeff=(vNaught)*dcdx1*sqrt(cofx2)/(2*sqrt(cofx1));
x1coeff=(vNaught^2)*dxstar1dx*dcdx2*sqrt(cofx1)/(2*sqrt(cofx2));
x2coeff=(vNaught^2)*dxstar2dx*dcdx1*sqrt(cofx2)/(2*sqrt(cofx1));
Emptycoeff=(vNaught^2)*dcdx1*sqrt(cofx2)/(2*sqrt(cofx1))*dcdx2*sqrt(cofx1)/(2*sqrt(cofx2));


X1T2t1t2coeff=(vNaught)*dt2dx*sqrt(cofx1*cofx2);
T1X2t1t2coeff=(vNaught)*dt1dx*sqrt(cofx1*cofx2);

X1X2t1coeff=(vNaught^2)*dt1dx*sqrt(cofx1)/(2*sqrt(cofx2));
X1X2t2coeff=(vNaught^2)*dt2dx*sqrt(cofx2)/(2*sqrt(cofx1));
X1X2t1x2coeff=(vNaught^2)*dt1dx*dxstar2dx*sqrt(cofx1*cofx2);
X1X2x1t2coeff=(vNaught^2)*dt2dx*dxstar1dx*sqrt(cofx1*cofx2);
X1X2t1t2coeff=(vNaught^2)*dt2dx*dt1dx*sqrt(cofx1*cofx2);




                AAAAA(m, n)=cofxsquaredterm*(Emptycoeff*real(ILower1(end))+real(IttLower1(end))+x1x2coeff*real(IxxLower1(end))+x1t2coeff*real(IxtLower1(end))+t1x2coeff*real(ItxLower1(end))+x1coeff*real(Ix1Lower1(end))+x2coeff*real(Ix2Lower1(end))+t1coeff*real(It1Lower1(end))+t2coeff*real(It2Lower1(end)));
               
                AAAAALR(m, n)=cofxsquaredterm*(Emptycoeff*real(ILRLower1(end))+real(ILRttLower1(end))+x1x2coeff*real(ILRxxLower1(end))+x1t2coeff*real(ILRxtLower1(end))+t1x2coeff*real(ILRtxLower1(end))+x1coeff*real(ILRx1Lower1(end))+x2coeff*real(ILRx2Lower1(end))+t1coeff*real(ILRt1Lower1(end))+t2coeff*real(ILRt2Lower1(end)));

      
                AAAAA2(m, n)=cofxsquaredterm*(X1T2t1t2coeff*real(IttLower1(end))+T1X2t1t2coeff*real(IttLower1(end))...
                    +X1X2t1coeff*real(It1Lower1(end))+X1X2t2coeff*real(It2Lower1(end))...
                    +X1X2t1x2coeff*real(ItxLower1(end))+X1X2x1t2coeff*real(IxtLower1(end))...
                    +X1X2t1t2coeff*real(IttLower1(end)));
                
                AAAAALR2(m, n)=cofxsquaredterm*(X1T2t1t2coeff*real(ILRttLower1(end))+T1X2t1t2coeff*real(ILRttLower1(end))...
                    +X1X2t1coeff*real(ILRt1Lower1(end))+X1X2t2coeff*real(ILRt2Lower1(end))...
                    +X1X2t1x2coeff*real(ILRtxLower1(end))+X1X2x1t2coeff*real(ILRxtLower1(end))...
                    +X1X2t1t2coeff*real(ILRttLower1(end)));
        end
    end
    
    
    dlmwrite([savefilepath,num2str(Region1PH),'/',num2str(Region2PH),'/',num2str(MassPH,'%10.5e\n'),TypePH,'AAAAAUpper.txt'],AAAAA,'precision','%30.25f','delimiter','\t')
    dlmwrite([savefilepath,num2str(Region1PH),'/',num2str(Region2PH),'/',num2str(MassPH,'%10.5e\n'),TypePH,'AAAAALRUpper.txt'],AAAAALR,'precision','%30.25f','delimiter','\t')
    
    
       dlmwrite([savefilepath,num2str(Region1PH),'/',num2str(Region2PH),'/',num2str(MassPH,'%10.5e\n'),TypePH,'AAAAA2Upper.txt'],AAAAA2,'precision','%30.25f','delimiter','\t')
    dlmwrite([savefilepath,num2str(Region1PH),'/',num2str(Region2PH),'/',num2str(MassPH,'%10.5e\n'),TypePH,'AAAAALR2Upper.txt'],AAAAALR2,'precision','%30.25f','delimiter','\t')
    
    
end
end
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    