// Test of iso226.sci
//
// ---------------------------------
// Federico Miyara
// 2025-07-03

clear, close(winsid())

// Load function
exec iso226.sci;

// Example 1
Ln1 = [20 40 60 80];
f =  [63 125 250 500 1000 2000 4000 8000];
v = 2023;
di = 1;
Lp1 = iso226(Ln1, f, v, di)


// Example 2
Lp2 = [20 40 60 80];
f =  [63 125 250 500 1000 2000 4000 8000];
v = 2023;
di = -1;
Ln2 = iso226(Lp2, f, v, di)

// ----------------------------
// ISO R226 : 1961
// ----------------------------

// Example 3
Ln3 = (1:10)*12;
f = 20*10^(0:0.05:log10(15000/20)); 
v = 1961;
di = 1;
[Lp3, Lt3]  = iso226(Ln3, f, v, di);

// Example 4
Lp4 = (1:10)*12;
f = 20*10^(0:0.05:log10(15000/20)); 
v = 1961;
di = -1;
Ln4 = iso226(Lp4, f, v, di);

scf(1);
clf(1);
gcf().figure_position = [200,80];
gcf().figure_size = [475,650];
subplot(2,1,1)
// Lp, parameter Ln
plot("ln", f, Lp3)
plot("ln", f, Lt3, "--r")
title("Isophone contours according to ISO R226:1961")
xlabel("f  (Hz)")
ylabel("Lp  (dB)")
gca().data_bounds = [20 -10; 20000 130];
gca().tight_limits = ["on","on","on"];
gca().grid = [1 1];
subplot(2,1,2)
// Ln, parpametro Lp
plot("ln", f, Ln4)
title("Equal sound pressure level contours according to ISO 226:1961")
xlabel("f  (Hz)")
ylabel("Ln  (fon)")
gca().data_bounds = [20 0; 20000 130];
gca().tight_limits = ["on","on","on"];
gca().grid = [1 1];


// ----------------------------
// ISO 226 : 1987
// ----------------------------

// Example 5
Ln5 = (1:10)*12;
f = 20*10^(0:0.05:log10(15000/20)); 
v = 1987;
di = 1;
[Lp5, Lt5]  = iso226(Ln5, f, v, di);

// Example 6
Lp6 = (1:10)*12;
f = 20*10^(0:0.05:log10(15000/20)); 
v = 1987;
di = -1;
Ln6 = iso226(Lp6, f, v, di);

scf(2);
clf(2);
gcf().figure_position = [220,60];
gcf().figure_size = [475,650];
subplot(2,1,1)
// Lp, parameter Ln
plot("ln", f, Lp5)
plot("ln", f, Lt5, "--r")
title("Isophone contours according to ISO 226:1987")
xlabel("f  (Hz)")
ylabel("Lp  (dB)")
gca().data_bounds = [20 -10; 20000 130];
gca().tight_limits = ["on","on","on"];
gca().grid = [1 1];
subplot(2,1,2)
// Ln, parameter Lp
plot("ln", f, Ln6)
title("Equal sound pressure level contours according to ISO 226:1987")
xlabel("f  (Hz)")
ylabel("Ln  (fon)")
gca().data_bounds = [20 0; 20000 130];
gca().tight_limits = ["on","on","on"];
gca().grid = [1 1];

// ----------------------------
// ISO 226 : 2003
// ----------------------------

// Example 7
Ln7 = (1:10)*12;
f = 20*10^(0:0.05:log10(15000/20)); 
v = 2003;
di = 1;
[Lp7, Lt7]  = iso226(Ln7, f, v, di);

// Example 8
Lp8 = (1:10)*12;
f = 20*10^(0:0.05:log10(15000/20)); 
v = 2003;
di = -1;
Ln8 = iso226(Lp8, f, v, di);

scf(3);
clf(3);
gcf().figure_position = [240,40];
gcf().figure_size = [475,650];
subplot(2,1,1)
// Lp, parameter Ln
plot("ln", f, Lp7)
plot("ln", f, Lt7, "--r")
title("Isophone contours according to ISO 226:2003")
xlabel("f  (Hz)")
ylabel("Lp  (dB)")
gca().data_bounds = [20 -10; 20000 130];
gca().tight_limits = ["on","on","on"];
gca().grid = [1 1];
subplot(2,1,2)
// Ln, parameter Lp
plot("ln", f, Ln8)
title("Equal sound pressure level contours according to ISO 226:2003")
xlabel("f  (Hz)")
ylabel("Ln  (fon)")
gca().data_bounds = [20 0; 20000 130];
gca().tight_limits = ["on","on","on"];
gca().grid = [1 1];

// ----------------------------
// ISO 226 : 2023
// ----------------------------

// Example 9
Ln9 = (1:10)*12;
f = 20*10^(0:0.05:log10(15000/20)); 
v = 2023;
di = 1;
[Lp9, Lt9]  = iso226(Ln9, f, v, di);

// Example 10
Lp10 = (1:10)*12;
f = 20*10^(0:0.05:log10(15000/20)); 
v = 2023;
di = -1;
Ln10 = iso226(Lp10, f, v, di);

scf(4);
clf(4);
gcf().figure_position = [260,20];
gcf().figure_size = [475,650];
subplot(2,1,1)
// Lp, parameter Ln
plot("ln", f, Lp9)
plot("ln", f, Lt9, "--r")
title("Isophone contours according to ISO 226:2023")
xlabel("f  (Hz)")
ylabel("Lp  (dB)")
gca().data_bounds = [20 -10; 20000 130];
gca().tight_limits = ["on","on","on"];
gca().grid = [1 1];
subplot(2,1,2)
// Ln, parameter Lp
plot("ln", f, Ln10)
title("Equal sound pressure level contours according to ISO 226:2023")
xlabel("f  (Hz)")
ylabel("Ln  (fon)")
gca().data_bounds = [20 0; 20000 130];
gca().tight_limits = ["on","on","on"];
gca().grid = [1 1];
