(* This computes the modes in chifunctions5.m *)

ClearAll

SetDirectory["/Users/zachscofield/OneDrive - Wake Forest University/Research/Simpson's"]

xb = 1/10;

x0=(101*xb)/100

m=1;

t=0;

xc = ((-xb)/2)+Sqrt[1-(3*xb^2)/4];

f = -(1/x)*(x-xc)*(x-xb)*(x+xc+xb);

fp = D[f,x];

κb=(1/(2*xb))*(xc-xb)*(xc+(2*xb))

κc=(1/(2*xc))*(xc-xb)*((2*xc)+xb)

κN=(1/(2(xc+xb)))*((2*xc)+xb)*(xc+(2*xb))

xstar = (1/(2*κb))*Log[(x-xb)/(xc-xb)]-(1/(2*κc))*Log[-(x-xc)/(xc-xb)]+(1/(2*κN))*Log[(x+xc+xb)/(xc+(2*xb))]-(xc/(4*xb*κb))*Log[((2*xc)+xb)/(xc+(2*xb))]-(xb*xc)/(2*(xc-xb))*Log[(xb/xc)];

x1a=(((10^9)-1)*xc)/(10^9)

ct1a=Exp[I*omegap*xstar]/.{x-> x1a};

cdt1a=((1/f) *I*omegap*Exp[I*omegap*xstar])/.{x-> x1a};

chiR1 = Table[{NDSolve[{chi''[x]+ ((1/f)*fp*chi'[x])+((i^2)/(f^2)-(m^2)/(f))*chi[x]==0, chi[x1a] == ct1a/.{omegap-> i}, chi'[x1a]==cdt1a/.{omegap-> i}}, chi, {x, x1a,x0}]},{i,1/10000,1/100,1/100000}];

chi[r]/.Part[chiR1[[100]],1]

chiR2 = Table[{NDSolve[{chi''[x]+ ((1/f)*fp*chi'[x])+((i^2)/(f^2)-(m^2)/(f))*chi[x]==0, chi[x1a] == ct1a/.{omegap-> i}, chi'[x1a]==cdt1a/.{omegap-> i}}, chi, {x, x1a,x0}]},{i,1/100,1,1/1000}];

x1b=(9999*xc)/10000

ct1b=Exp[I*omegap*xstar]/.{x-> x1b};

cdt1b=((1/f) *I*omegap*Exp[I*omegap*xstar])/.{x-> x1b};

chiR3 = Table[{NDSolve[{chi''[x]+ ((1/f)*fp*chi'[x])+((i^2)/(f^2)-(m^2)/(f))*chi[x]==0, chi[x1b] == ct1b/.{omegap-> i}, chi'[x1b]==cdt1b/.{omegap-> i}}, chi, {x, x1b,x0}]},{i,1,6,5/1000}];

chiR4 = Table[{NDSolve[{chi''[x]+ ((1/f)*fp*chi'[x])+((i^2)/(f^2)-(m^2)/(f))*chi[x]==0, chi[x1b] == ct1b/.{omegap-> i}, chi'[x1b]==cdt1b/.{omegap-> i}}, chi, {x, x1b,x0}]},{i,6,5000/100,4/100}];

chiR5 = Table[{NDSolve[{chi''[x]+ ((1/f)*fp*chi'[x])+((i^2)/(f^2)-(m^2)/(f))*chi[x]==0, chi[x1b] == ct1b/.{omegap-> i}, chi'[x1b]==cdt1b/.{omegap-> i}}, chi, {x, x1b,x0}]},{i,5000/100,10000/100,5/100}];

ct2=Exp[-I*omegap*xstar]/.{x-> x1a};

cdt2=((-1/f) *I*omegap*Exp[-I*omegap*xstar])/.{x-> x1a};

chiL1 = Table[{NDSolve[{chi''[x]+ ((1/f)*fp*chi'[x])+((i^2)/(f^2)-(m^2)/(f))*chi[x]==0, chi[x1a] == ct2/.{omegap-> i}, chi'[x1a]==cdt2/.{omegap-> i}}, chi, {x, x1a,x0}]},{i,1/10000,1/100,1/100000}];

chiL2 = Table[{NDSolve[{chi''[x]+ ((1/f)*fp*chi'[x])+((i^2)/(f^2)-(m^2)/(f))*chi[x]==0, chi[x1a] == ct2/.{omegap-> i}, chi'[x1a]==cdt2/.{omegap-> i}}, chi, {x, x1a,x0}]},{i,1/100,1,1/1000}];

ct2b=Exp[-I*omegap*xstar]/.{x-> x1b};

cdt2b=((-1/f) *I*omegap*Exp[-I*omegap*xstar])/.{x-> x1b};

chiL3 = Table[{NDSolve[{chi''[x]+ ((1/f)*fp*chi'[x])+((i^2)/(f^2)-(m^2)/(f))*chi[x]==0, chi[x1b] == ct2b/.{omegap-> i}, chi'[x1b]==cdt2b/.{omegap-> i}}, chi, {x, x1b,x0}]},{i,1,6,5/1000}];

chiL4 = Table[{NDSolve[{chi''[x]+ ((1/f)*fp*chi'[x])+((i^2)/(f^2)-(m^2)/(f))*chi[x]==0, chi[x1b] == ct2b/.{omegap-> i}, chi'[x1b]==cdt2b/.{omegap-> i}}, chi, {x, x1b,x0}]},{i,6,5000/100,4/100}];

chiL5 = Table[{NDSolve[{chi''[x]+ ((1/f)*fp*chi'[x])+((i^2)/(f^2)-(m^2)/(f))*chi[x]==0, chi[x1b] == ct2b/.{omegap-> i}, chi'[x1b]==cdt2b/.{omegap-> i}}, chi, {x, x1b,x0}]},{i,5000/100,10000/100, 5/100}];

(* Save["chifunctions5.m",{chiR1,chiR2,chiR3,chiR4,chiR5, chiL1,chiL2,chiL3,chiL4,chiL5}] *)