In[1]:= ClearAll[f];
refrWater = InterpolatingFunction[{{182.,200000.}},{5,7,0,{455},{4},0,0,0,0,Automatic,{},{},False},{{182.,191.5,200.9,210.4,219.9,229.4,238.8,248.29999999999998,257.79999999999995,267.2,276.7,286.2,295.59999999999997,305.09999999999997,314.59999999999997,324.1,333.5,343.,352.5,361.9,371.40000000000003,380.90000000000003,390.29999999999995,399.8,409.3,418.8,428.20000000000005,437.7,447.2,456.6,466.1,475.6,485.,494.5,504.,513.5,522.9,532.4,541.9000000000001,551.3000000000001,560.8,570.3000000000001,579.7,589.1999999999999,598.7,608.1999999999999,617.6,627.1,636.6,646.,655.5,665.,674.4,683.9,693.4,702.9,712.3000000000001,721.8,731.3,740.7,750.1999999999999,759.7,769.1,778.5999999999999,788.1,797.6,807.,816.5,826.,835.4,844.9,854.4000000000001,863.8000000000001,873.3,882.8000000000001,892.3,901.6999999999999,911.2,920.6999999999999,930.1,939.6,949.1,958.5,968.,977.5,987.,996.4,1006.,1014.9999999999999,1025.,1034.,1044.,1053.,1063.,1072.,1082.,1091.,1101.,1110.,1120.,1129.,1200.,1400.,1600.,1800.,2000.,2020.,2041.,2062.,2083.,2105.,2128.,2151.,2174.,2198.,2222.,2247.,2273.,2299.,2326.,2353.,2381.,2410.,2439.,2469.,2500.,2506.,2513.,2519.,2525.,2532.,2538.,2545.,2551.,2558.,2564.,2571.,2577.,2584.,2591.,2597.,2604.,2611.,2618.,2625.,2632.,2639.,2646.,2653.,2660.,2667.,2674.,2681.,2688.,2695.,2703.,2710.,2717.,2725.,2732.,2740.,2747.,2755.,2762.,2770.,2778.,2786.,2793.,2801.,2809.,2817.,2825.,2833.,2841.,2849.,2857.,2865.,2874.,2882.,2890.,2899.,2907.,2915.,2924.,2933.,2941.,2950.,2959.,2967.,2976.,2985.,2994.,3003.,3012.,3021.,3030.,3040.,3049.,3058.,3067.,3077.,3086.,3096.,3106.,3115.,3125.,3135.,3145.,3155.,3165.,3175.,3185.,3195.,3205.,3215.,3226.,3236.,3247.,3257.,3268.,3279.,3289.,3300.,3311.,3322.,3333.,3344.,3356.,3367.,3378.,3390.,3401.,3413.,3425.,3436.,3448.,3460.,3472.,3484.,3497.,3509.,3521.,3534.,3546.,3559.,3571.,3584.,3597.,3610.,3623.,3636.,3650.,3663.,3676.,3690.,3704.,3717.,3731.,3745.,3759.,3774.,3788.,3802.,3817.,3831.,3846.,3861.,3876.,3891.,3906.,3922.,3937.,3953.,3968.,3984.,4000.,4016.,4032.,4049.,4065.,4082.,4098.,4115.,4132.,4149.,4167.,4184.,4202.,4219.,4237.,4255.,4274.,4292.,4310.,4329.,4348.,4367.,4386.,4405.,4425.,4444.,4464.,4484.,4505.,4525.,4545.,4566.,4587.,4608.,4630.,4651.,4673.,4695.,4717.,4739.,4762.,4785.,4808.,4831.,4854.,4878.,4902.,4926.,4950.,4975.,5000.,5025.,5051.,5076.,5102.,5128.,5155.,5181.,5208.,5236.,5263.,5291.,5319.,5348.,5376.,5405.,5435.,5464.,5495.,5525.,5556.,5587.,5618.,5650.,5682.,5714.,5747.,5780.,5814.,5848.,5882.,5917.,5952.,5988.,6024.,6061.,6098.,6135.,6173.,6211.,6250.,6289.,6329.,6369.,6410.,6452.,6494.,6536.,6579.,6623.,6667.,6700.,6800.,6900.,7000.,7100.,7200.,7300.,7400.,7500.,7600.,7700.,7800.,7900.,8000.,8200.,8400.,8600.,8800.,9000.,9200.,9400.,9600.,9800.,10000.,10500.,11000.,11500.,12000.,12500.,13000.,13500.,14000.,14500.,15000.,15500.,16000.,16500.,17000.,17500.,18000.,18500.,19000.,19500.,20000.,21000.,22000.,23000.,24000.,25000.,26000.,27000.,28000.,29000.,30000.,32000.,34000.,36000.,38000.,40000.,42000.,44000.,46000.,48000.,50000.,60000.,70000.,80000.,90000.,100000.,110000.,120000.,130000.,140000.,150000.,160000.,170000.,180000.,190000.,200000.}},{Developer`PackedArrayForm,{0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273,274,275,276,277,278,279,280,281,282,283,284,285,286,287,288,289,290,291,292,293,294,295,296,297,298,299,300,301,302,303,304,305,306,307,308,309,310,311,312,313,314,315,316,317,318,319,320,321,322,323,324,325,326,327,328,329,330,331,332,333,334,335,336,337,338,339,340,341,342,343,344,345,346,347,348,349,350,351,352,353,354,355,356,357,358,359,360,361,362,363,364,365,366,367,368,369,370,371,372,373,374,375,376,377,378,379,380,381,382,383,384,385,386,387,388,389,390,391,392,393,394,395,396,397,398,399,400,401,402,403,404,405,406,407,408,409,410,411,412,413,414,415,416,417,418,419,420,421,422,423,424,425,426,427,428,429,430,431,432,433,434,435,436,437,438,439,440,441,442,443,444,445,446,447,448,449,450,451,452,453,454,455},{1.4677097441493,1.4404331506401,1.4223657629504,1.4089720452599,1.3986433047197,1.390400109458,1.383721037604,1.37808831852,1.3733214421885,1.3692768869158,1.3657352327502,1.3626415791521,1.3599451690455,1.3575283523326,1.3553733211095,1.3534412070092,1.3517178003209,1.3501406237405,1.3487073377948,1.3474130030283,1.346214937405,1.345114563197,1.3441110456812,1.3431735819155,1.342305043532,1.3414982764576,1.3407546807519,1.3400530171597,1.3393966700296,1.3387876724556,1.338209340659,1.3376650208921,1.3371569860573,1.3366717871842,1.3362125993297,1.3357772766023,1.3353681293216,1.3349747209786,1.3345999349096,1.3342460158649,1.3339041405322,1.3335769734208,1.3332666628587,1.3329656231228,1.3326763163567,1.3323979087118,1.332132407632,1.3318734648367,1.3316233062778,1.331383837792,1.3311494136837,1.33092210909,1.3307037367741,1.3304892114579,1.330280478546,1.3300771547067,1.3298809446661,1.3296873432698,1.3294981523223,1.3293150081137,1.3291337434839,1.3289560721242,1.3287835722159,1.3286123499273,1.3284440453935,1.3282784588304,1.3281171057973,1.3279563771328,1.3277978345172,1.3276429567775,1.3274882997613,1.3273353791879,1.3271856450395,1.327035784463,1.3268872767838,1.3267400061551,1.3265953951297,1.3264502626387,1.326306052485,1.3261641741362,1.3260215180451,1.3258795099981,1.325739552157,1.3255985867561,1.325458028092,1.325317801732,1.3251793082663,1.3250380620977,1.3249057611921,1.3247588255499,1.324626584498,1.3244795867846,1.3243471777947,1.3241998728782,1.3240670820459,1.3239192392189,1.3237858647693,1.3236372659102,1.3235031166389,1.3233535545788,1.3232184484642,1.324,1.321,1.317,1.312,1.303,1.301,1.301,1.3,1.298,1.298,1.296,1.295,1.294,1.293,1.291,1.289,1.287,1.285,1.282,1.28,1.277,1.274,1.27,1.265,1.261,1.26,1.259,1.257,1.256,1.255,1.254,1.252,1.25,1.249,1.247,1.246,1.243,1.241,1.24,1.238,1.235,1.232,1.23,1.227,1.224,1.221,1.218,1.214,1.21,1.205,1.2,1.195,1.191,1.185,1.179,1.172,1.166,1.157,1.149,1.144,1.139,1.138,1.138,1.139,1.141,1.144,1.149,1.154,1.158,1.161,1.165,1.171,1.177,1.183,1.191,1.199,1.212,1.22,1.233,1.246,1.258,1.271,1.282,1.293,1.305,1.317,1.329,1.342,1.353,1.364,1.376,1.386,1.398,1.407,1.417,1.426,1.434,1.442,1.45,1.457,1.465,1.471,1.476,1.48,1.483,1.486,1.487,1.487,1.487,1.486,1.485,1.482,1.479,1.477,1.474,1.472,1.467,1.464,1.461,1.457,1.454,1.451,1.448,1.444,1.441,1.437,1.434,1.431,1.427,1.425,1.421,1.418,1.415,1.413,1.41,1.407,1.405,1.403,1.4,1.398,1.396,1.394,1.392,1.39,1.388,1.387,1.385,1.383,1.382,1.379,1.378,1.377,1.375,1.374,1.372,1.371,1.37,1.369,1.367,1.366,1.365,1.363,1.361,1.361,1.36,1.358,1.358,1.357,1.355,1.354,1.353,1.352,1.351,1.35,1.349,1.348,1.348,1.347,1.346,1.345,1.344,1.344,1.343,1.342,1.341,1.34,1.34,1.338,1.337,1.337,1.335,1.334,1.334,1.333,1.332,1.332,1.331,1.33,1.33,1.33,1.329,1.329,1.329,1.328,1.328,1.327,1.327,1.327,1.327,1.327,1.326,1.326,1.326,1.325,1.325,1.325,1.325,1.325,1.325,1.324,1.324,1.323,1.322,1.322,1.321,1.32,1.319,1.318,1.318,1.317,1.316,1.314,1.313,1.311,1.31,1.308,1.306,1.304,1.302,1.299,1.297,1.294,1.291,1.288,1.285,1.282,1.278,1.275,1.271,1.267,1.262,1.256,1.251,1.247,1.242,1.241,1.241,1.247,1.265,1.289,1.311,1.332,1.349,1.354,1.356,1.354,1.35,1.345,1.341,1.337,1.333,1.33,1.326,1.324,1.322,1.329,1.324,1.321,1.317,1.314,1.312,1.309,1.307,1.304,1.302,1.299,1.297,1.294,1.291,1.286,1.281,1.275,1.269,1.262,1.255,1.247,1.239,1.229,1.218,1.185,1.153,1.126,1.111,1.123,1.146,1.177,1.21,1.241,1.27,1.297,1.325,1.351,1.376,1.401,1.423,1.443,1.461,1.476,1.48,1.487,1.5,1.511,1.521,1.531,1.539,1.545,1.549,1.551,1.551,1.546,1.536,1.527,1.522,1.519,1.522,1.53,1.541,1.555,1.587,1.703,1.821,1.886,1.924,1.957,1.966,2.004,2.036,2.056,2.069,2.081,2.094,2.107,2.119,2.13}},{Automatic}];

(* fundamental constants *)
eR = 78; (*relative permittivity, bulk water*)
e0 = 8.854*^-12 ;(* vacuum permittivity F/m *)
kB = 1.3806*^-23;(*boltzmann constant, J/K*)
elec = 1.6021*^-19;(*elementary charge electron, C*)
T = 298; (*temp Kelvin*)
fara=96488;(*Faraday constant, C/mol*)
avo=6.023*^23; (*Avogadro's number, /mol*)
c = 2.998*^8 ;(* speed of light, m/s *)

(*********** SET: OPTICAL PARAMETERS *************)
wl1 = 532; (* visible wavelength, nm *)


(* wavenumbers of broadband beam *)
wLow = 2700; 
wHigh = 3700;
wRes = .5;
w = Range[wLow,wHigh,wRes];

(* Incident angles in degrees *)
theta1 = 60; (* visible *)
theta2 = 56; (* IR *)

(*calculate: dKZ for specific wavenumbers given above parameters*)
dKZfromWN[wn_]:=( 2*Pi*c/(wl1*10^-9))/c*Sqrt[refrWater[wl1]^2 - Sin[theta1 °]^2]+ ( 2*Pi*c/((10^7/wn)*10^-9))/c*Sqrt[refrWater[10^7/wn]^2 - Sin[theta2 °]^2] + (2*Pi*c/((10^7/((10^7/wl1)+wn))*10^-9))/c*Sqrt[refrWater[(10^7/((10^7/wl1)+wn))]^2 - Sin[Pi/2-ArcTan[((2*Pi/(wl1*10^-9))*Cos[theta1 °]+ (2*Pi/((10^7/wn)*10^-9))*Cos[theta2 °])/((2*Pi/(wl1*10^-9))*Sin[theta1 °]+ (2*Pi/((10^7/wn)*10^-9))*Sin[theta2 °])]]^2];

(********** SET: peak parameters ************)

(* PEAK 1 OH 3200 and 3250 *)
(* chi2 *)
w1 = 3200; (* position of peak 1 *)
chi2w1 = 10; (* magnitude of peak 1 chi2 *)
phiChi2w1= Pi; (*phase of peak 1 chi 2 for DODA*)
phiChi2w1CNC= 0; (*phase of peak 1 chi 2 for DODA/CNC*)
gammaChi2w1 = 120; (*damping of peak 1 chi 2 *)
(* chi3 *)
w1chi3 = 3250; (* shifted position of peak 1 chi3 *)
chi3w1 = 500; (* magnitude of peak 1 chi3 *)
phiChi3w1 = 87/180*Pi; (* phase of peak 1 chi3 for DODA*)
phiChi3w1CNC = 77/180*Pi; (* phase of peak 1 chi3 for DODA/CNC*)
gammaChi3w1 = 150; (*damping of peak 1 chi 2 *)

(* PEAK2 OH 3400 and 3450 *)
(* chi2 *)
w2 = 3400; (* position of peak 2 *)
chi2w2 = 10;(* magnitude of peak 2 chi2 *)
phiChi2w2= Pi; (* phase of peak 2 chi2 for DODA*)
phiChi2w2CNC= 0; (* phase of peak 2 chi2 for DODA/CNC*)
gammaChi2w2 = 120; (*damping of peak 1 chi 2 *)
(* chi3 *)
w2chi3 = 3450; (* shifted position of peak 2 chi3 *)
chi3w2 =150; (* magnitude of peak 2 chi3 *)
phiChi3w2 = 87/180*Pi; (* phase of peak 2 chi3 for DODA *)
phiChi3w2CNC = 77/180*Pi; (* phase of peak 2 chi3 for DODA/CNC *)
gammaChi3w2 = 150; (*damping of peak 1 chi 2 *)

(* PEAK 3 DODA Gauche defects*)
(* chi2 *)
w3 = 2854; (* position of peak 3 *)
chi2w3 = 1.2;(* magnitude of peak 3 chi2 *)
phiChi2w3= Pi; (* phase of peak 3 chi2 *)
gammaChi2w3 = 10.9; (*damping of peak 3 chi 2 *)

(* PEAK 4 DODA CH3 (ss) *)
(* chi2 *)
w4 = 2875; (* position of peak 4 *)
chi2w4 = 4;(* magnitude of peak 4 chi2 *)
phiChi2w4= Pi; (* phase of peak 4 chi2 *)
gammaChi2w4 = 8.8; (*damping of peak 4 chi 2 *)

(* PEAK 5 DODA a-CH2 (ss) *)
(* chi2 *)
w5 = 2887; (* position of peak 5 *)
chi2w5 = 1.4;(* magnitude of peak 5 chi2 *)
phiChi2w5= 0; (* phase of peak 5 chi2 *)
gammaChi2w5 = 7.9; (*damping of peak 5 chi 2 *)

(* PEAK 6 DODA a-CH2 (as)*)
(* chi2 *)
w6 = 2911; (* position of peak 6 *)
chi2w6 = 1.4;(* magnitude of peak 6 chi2 *)
phiChi2w6= Pi; (* phase of peak 6 chi2 *)
gammaChi2w6 = 13.1; (*damping of peak 6 chi 2 *)

(* PEAK 7 DODA CH3 (FR)*)
(* chi2 *)
w7 = 2934; (* position of peak 7 *)
chi2w7 = 3.1;(* magnitude of peak 7 chi2 *)
phiChi2w7= Pi; (* phase of peak 7 chi2 *)
gammaChi2w7 = 8.7; (*damping of peak 7 chi 2 *)

(* PEAK 8 DODA CH3 (as)*)
(* chi2 *)
w8 = 2956; (* position of peak 8 *)
chi2w8 = 3.2;(* magnitude of peak 8 chi2 *)
phiChi2w8= 0; (* phase of peak 8 chi2 *)
gammaChi2w8 = 10.9; (*damping of peak 8 chi 2 *)

(* PEAK 9 DODA CH3 (as)*)
(* chi2 *)
w9 = 2992; (* position of peak 9 *)
chi2w9 = 0.2;(* magnitude of peak 9 chi2 *)
phiChi2w9= 0; (* phase of peak 9 chi2 *)
gammaChi2w9 = 9.3; (*damping of peak 9 chi 2 *)

(*********** SECTION 1: ionic strength constant ********************)

(************ SET Ionic strength [M] for DODA ****************)
conc =2.61*10^(-7); 
(************ SET Ionic strength [M] for DODA/CNC ****************)
concCNC =3.98*10^(-6); 

(********** SET: surface charge density [C/m^2] for DODA******)
sigmas = {0.04,0.001};
(********** SET: surface charge density [C/m^2] for DODA/CNC******)
sigmasCNC = {-0.04,-0.001};

(* calculate: Debye Length [m] *)
debye = Sqrt[(eR*e0*8.314*T)/(2*(fara^2)*conc*1000)];
kappa = 1/(1*debye);

(* calculate: Debye Length [m] for DODA/CNC*)
debyeCNC = Sqrt[(eR*e0*8.314*T)/(2*(fara^2)*concCNC*1000)];
kappaCNC = 1/(1*debyeCNC);

(* calculate: potential by Gouy Chapman model for DODA*)
potentials = ((2*kB*T)/elec)*ArcSinh[sigmas/Sqrt[8000*kB*T*avo*eR*e0*conc]];
Print["potentials DODA= ", potentials*1000, " mV"];
Print["Debye length DODA= ",debye*1*^9, " nm"];

(* calculate: potential by Gouy Chapman model for DODA/CNC*)
potentialsCNC = ((2*kB*T)/elec)*ArcSinh[sigmasCNC/Sqrt[8000*kB*T*avo*eR*e0*concCNC]];
Print["potentials DODA/CNC = ", potentialsCNC*1000, " mV"];
Print["Debye length DODA/CNC = ",debyeCNC*1*^9, " nm"];
Print["dKZ = ",dKZfromWN[3250]*1*^-9, "nm-1"];
Print["kappa DODA = ",1/debye*1*^-9, " nm-1"];
Print["kappa DODA/CNC = ",1/debyeCNC*1*^-9, " nm-1"];
Print["phichi3w1,2 DODA = ",ArcTan[kappa,dKZfromWN[3250]]/Pi*180, " °"];
Print["phichi3w1,2 DODA/CNC = ",ArcTan[kappaCNC,dKZfromWN[3250]]/Pi*180, " °"];

(**********PLOT 1: DDOA case ************)
fONeb[pot_]:= (chi2w1*Exp[I*phiChi2w1])/(w-w1+gammaChi2w1*I)+(chi2w2*Exp[I*phiChi2w2])/(w-w2+gammaChi2w2*I)+(chi2w3*Exp[I*phiChi2w3])/(w-w3+gammaChi2w3*I)+(chi2w4*Exp[I*phiChi2w4])/(w-w4+gammaChi2w4*I)+(chi2w5*Exp[I*phiChi2w5])/(w-w5+gammaChi2w5*I)+(chi2w6*Exp[I*phiChi2w6])/(w-w6+gammaChi2w6*I)+(chi2w7*Exp[I*phiChi2w7])/(w-w7+gammaChi2w7*I)+(chi2w8*Exp[I*phiChi2w8])/(w-w8+gammaChi2w8*I)+(chi2w9*Exp[I*phiChi2w9])/(w-w9+gammaChi2w9*I)+((chi3w1*Exp[I*phiChi3w1])/(w-w1chi3+gammaChi3w1*I)+(chi3w2*Exp[I*phiChi3w2])/(w-w2chi3+gammaChi3w2*I))*(kappa/(Sqrt[kappa^2+dKZfromWN[w]^2]))Exp[I*ArcTan[dKZfromWN[w]/kappa]]*pot;


ydata = Abs[Map[fONeb,potentials]]^2;
plotdata=Inner[List,w,#,List]&/@ ydata;
ListLinePlot[plotdata,Frame->True,BaseStyle->{FontFamily->"Arial",FontSize->16},ImageSize->{450,270},PlotRange->{{2800,3000},{-0.02,0.5}},PlotLegends-> sigmas,PlotLabel->"DODA at the air/water interface"]

(**********PLOT 2: DODA/CNC case ************)
fONebCNC[potCNC_]:= (chi2w1*Exp[I*phiChi2w1CNC])/(w-w1+gammaChi2w1*I)+(chi2w2*Exp[I*phiChi2w2CNC])/(w-w2+gammaChi2w2*I)+(chi2w3*Exp[I*phiChi2w3])/(w-w3+gammaChi2w3*I)+(chi2w4*Exp[I*phiChi2w4])/(w-w4+gammaChi2w4*I)+(chi2w5*Exp[I*phiChi2w5])/(w-w5+gammaChi2w5*I)+(chi2w6*Exp[I*phiChi2w6])/(w-w6+gammaChi2w6*I)+(chi2w7*Exp[I*phiChi2w7])/(w-w7+gammaChi2w7*I)+(chi2w8*Exp[I*phiChi2w8])/(w-w8+gammaChi2w8*I)+(chi2w9*Exp[I*phiChi2w9])/(w-w9+gammaChi2w9*I)+((chi3w1*Exp[I*phiChi3w1CNC])/(w-w1chi3+gammaChi3w1*I)+(chi3w2*Exp[I*phiChi3w2CNC])/(w-w2chi3+gammaChi3w2*I))*(kappaCNC/(Sqrt[kappaCNC^2+dKZfromWN[w]^2]))Exp[I*ArcTan[dKZfromWN[w]/kappaCNC]]*potCNC;


ydataCNC = Abs[Map[fONebCNC,potentialsCNC]]^2;
plotdataCNC=Inner[List,w,#,List]&/@ ydataCNC;
ListLinePlot[plotdataCNC,Frame->True,BaseStyle->{FontFamily->"Arial",FontSize->16},ImageSize->{450,270},PlotRange->{{2800,3000},{-0.02,0.5}},PlotLegends-> sigmasCNC,PlotLabel->"DODA/CNC at the air/water interface"]