function [Lp, Lt] = iso226(Ln, f, v, di)
   // Sound pressure level from loudness level or viceversa
   //
   // Usage:
   //       Lp = iso226(Ln, f, v, di)
   // where
   //       Ln: Loudness level in phon or sound
   //           pressure level in dB (if di = -1) 
   //       f:  Frequency in Hz
   //       v:  Version of the ISO226 Standard: 
   //           1961, 1987, 2003, 2023
   //       di: Direction of computation. If -1,
   //           Lp and Ln are toggled
   //       Lp: Sound pressure level in dB or 
   //           loudness level in phon (if di = -1)
   //       Lt: Hearing threshold at f in dB
   //
   // NOTE: Ln may be a scalar or a column vector
   //       f may be a scalar or a row vector
   //       In such case Lp is a matrix with as
   //       many columns as values of f and as many 
   //       rows as values of Ln.
   //
   // --------------------------------
   // Author: Federico Miyara
   // Date:   2025-07-02   

   rhs = argn(2)
   
   // Missing arguments
   if rhs<4
      di = 1;
   end
   if rhs<3
      v = 2023;
   end
   
   // If v is not a valid date  
   if and(v<>[1961 1987 2003 2023])
      mprintf("Version not valid, using 2023\n")
      v = 2023;   
   end
   
   select v
   case 1961
      // Parameters oorresponding to Robinson 
      // and Dadson contours 
      f61 = [     20     25     30     35    40    45    50    55 ..   
                  60     70     80     90   100   120   140   160 ..   
                 180    200    250    300   350   400   500   600 ..
                 700    800    900   1000  1200  1400  1600  1800 ..
                2000   2500   3000   3500  4000  4500  5000  6000 ..
                7000   8000   9000  10000 11000 12000 13000 14000 ..
               15000];
      af61 = [-217.2 -167.2 -135.7 -113.2 -96.3 -83.2 -73.0 -65.1 ..
               -58.7  -49.1  -42.5  -37.4 -33.5 -27.0 -22.7 -19.4 ..
               -16.9  -14.7  -10.8   -8.1  -6.1  -4.7  -3.0  -1.8 ..
                -1.0   -0.5   -0.2    0.0   0.5   1.1   1.7   2.5 ..
                 3.3    5.3    6.9    7.9   7.9   7.1   5.3  -0.5 ..
                -7.5  -13.3  -16.5  -16.8 -14.8 -12.7 -13.9 -22.7 ..
               -43.0]; 
      bf61 = [ 3.669  3.145  2.801  2.538 2.325 2.154 2.021 1.918 ..
               1.836  1.723  1.652  1.603 1.570 1.512 1.473 1.444 ..
               1.422  1.404  1.362  1.325 1.290 1.259 1.205 1.155 ..
               1.109  1.064  1.028  1.000 0.972 0.953 0.937 0.927 ..
               0.924  0.928  0.937  0.946 0.954 0.963 0.973 1.011 ..
               1.075  1.159  1.242  1.314 1.377 1.450 1.566 1.777 ..
               2.146];
      cf61 = [  -927   -782   -677   -589  -515  -454  -407  -370 ..
                -341   -305   -285   -274  -269  -260  -253  -248 ..
                -244   -242   -231   -218  -202  -185  -151  -117 ..
                 -84    -50    -22      0    23    47    71    90 ..
                 100    118    120    113    98    80    59    14 ..
                 -35    -89   -145   -203  -269  -350  -454  -591 ..
                 772]/1e5;
      Luf61 = [ 74.3   65.0   58.1   52.7  48.4  44.8  41.7  39.1 ..
                36.8   32.9   29.8   27.2  25.1  21.4  18.9  16.8 ..
                15.2   13.8   11.2    9.4   8.1   7.2   6.0   5.2 ..
                 4.7    4.4    4.3    4.2   3.9   3.4   2.6   1.8 ..
                 1.0   -1.2   -2.9   -3.9  -3.9  -3.0  -1.1   4.6 ..
                10.9   15.3   17.0   16.4  14.2  12.0  12.0  16.0 ..
                24.1];

      // Interpolation of the parameters af, bf, cf
      // at f
      aff = interp1(f61, af61, f, "spline");
      bff = interp1(f61, bf61, f, "spline");
      cff = interp1(f61, cf61, f, "spline");
      Lt =  interp1(f61, Luf61, f, "spline");

      if di==-1
         // Loudness level calculation
         // NOTE: To make the syntax more versatile
         //       here Lp stands for loudness level
         //       and Ln for sound pressure level 
         for k=1:length(Ln)
            Lp(k,:) = aff + bff*Ln(k) + cff*Ln(k).^2; 
         end
      else
         // Sound pressure level calculation
         // The mathematical model gives the loudness 
         // level as a second degree polynomial of 
         // the the sound pressure level so to get 
         // the sound pressure level from the loudness 
         // level it is necessary to solve the second
         // degree equation.
         for k=1:length(Ln)
            A = cff;
            B = bff;
            C = aff - Ln(k);
            Lp(k,:) = (-B + sqrt(B.^2 - 4*A.*C))/2./A;
            // Correct the values where A=0 y where the 
            // quadratic formula is not applicable
            Lp(k,find(A==0)) = Ln(k);
         end      
      end 
   
   case 1987
      // Parameters corresponding to ISO 226:1987 
      // One third octave normalized frequencies
      f87  = [   20    25  31.5    40    50    63    80   100 ..
                125   160   200   250   315   400   500   630 ..
                800  1000  1250  1600  2000  2500  3150  4000 ..
               5000  6300  8000 10000 12500];
      af87 = [2.347 2.190 2.050 1.879 1.724 1.597 1.512 1.466 ..
              1.426 1.394 1.372 1.344 1.304 1.256 1.203 1.135 ..
              1.062 1.000 0.967 0.943 0.932 0.933 0.937 0.952 ..
              0.974 1.027 1.135 1.266 1.501];
      bf87 = [ 5.61  5.27  4.81  4.04  3.38  2.86  2.59  2.57 ..
               2.56  2.55  2.54  2.48  2.29  2.01  1.62  1.11 ..
               0.52  0.00 -0.39 -0.67 -0.92 -1.05 -1.04 -0.88 ..
              -0.55  0.00  0.89  2.11  4.88] * 0.001;
      Tf87 = [ 74.3  65.0  56.3  48.4  41.7  35.5  29.8  25.1 ..
               20.7  16.8  13.8  11.2   8.9   7.2   6.0   5.0 ..
                4.4   4.2   3.7   2.6   1.0  -1.2  -3.6  -3.9 ..
               -1.1   6.6  15.3  16.4  11.6]; 

      // With these parameters it is possible to 
      // compute the loudness level Ln from the
      // sound pressure level Lp 
      //
      //               af(Lp - Tf)
      // Ln = 4.2 + -----------------
      //             1 + bf(Lp - Tf)  
      //
      // or the sound pressure level Ln from the
      // loudness level Lp 
      //
      //                Ln - 4.2
      // Lp = Tf + -------------------
      //            af - bf(Ln - 4.2)

      // Interpolation of the parameters af, bf, Tf
      // at f
      aff = interp1(f87, af87, f, "spline");
      bff = interp1(f87, bf87, f, "spline");
      Tff = interp1(f87, Tf87, f, "spline");
      Lt = Tff;
      
      if di==-1
         // Loudness level
         for k=1:length(Ln)
            Lp(k, :) = 4.2 + aff.*(Ln(k) - Tff) ./ ..
                       (1 + bff.*(Ln(k) - Tff));
         end      
      else
         // Sound pressure level
         for k=1:length(Ln)
            Lp(k, :) = Tff + (Ln(k) - 4.2) ./ ..
                       (aff - bff*(Ln(k) - 4.2));
         end
      end
      
   case 2003
      // Parameters corresponding to the contours
      // of Suzuki-Takeshima (?) and ISO 226:2003 
      f03  = [      20    25  31.5    40    50    63    80   100 ..
                   125   160   200   250   315   400   500   630 ..
                   800  1000  1250  1600  2000  2500  3150  4000 ..
                  5000  6300  8000 10000 12500];
      alfaf03 = [0.532 0.506 0.480 0.455 0.432 0.409 0.387 0.367 ..
                 0.349 0.330 0.315 0.301 0.288 0.276 0.267 0.259 ..
                 0.253 0.250 0.246 0.244 0.243 0.243 0.243 0.242 ..
                 0.242 0.245 0.254 0.271 0.301];
      LUf03 = [  -31.6 -27.2 -23.0 -19.1 -15.9 -13.0 -10.3  -8.1 ..
                  -6.2  -4.5  -3.1  -2.0  -1.1  -0.4   0.0   0.3 ..
                   0.5   0.0  -2.7  -4.1  -1.0   1.7   2.5   1.2 ..
                  -2.1  -7.1 -11.2 -10.7  -3.1]; 
      Tf03 = [    78.5  68.7  59.5  51.1  44.0  37.5  31.5  26.5 ..
                  22.1  17.9  14.4  11.4   8.6   6.2   4.4   3.0 ..
                   2.2   2.4   3.5   1.7  -1.3  -4.2  -6.0  -5.4 .. 
                  -1.5   6.0  12.6  13.9  12.3];

      // With these parameters it is possible to
      // calculate the loudness level Ln from 
      // the sound presssure level Lp: 
      //
      // Ln = 40*log10(Bf) + 94
      // 
      // where
      //
      // Bf = (0.4 * 10.^((Lp + LU - 90)/10)).^alfaf - .. 
      //      (0.4 * 10.^((Tf + LU - 90)/10)).^alfaf + 0.005135
      //
      // or the sound pressure level Lp corresponding
      // to a louness level Ln:
      //
      // Lp = (10/alfaf)* log(Af) - LU + 94
      //
      // where 
      //
      // Af = 4.47e-3*(10.^(Ln/40) - 1.15) + 
      //     (0.4*10.^((TTf+LU-90)/10)).^alfaf
      
      // Interpolation of the parameters alfaf, 
      // LUf, Tf at f
      alfaff = interp1(f03, alfaf03, f, "spline");
      LUff = interp1(f03, LUf03, f, "spline");
      Tff = interp1(f03, Tf03, f, "spline");
      Lt = Tff;
      
      if di==-1
         // Loudness level
         for k=1:length(Ln)
            Bff = (0.4 * 10.^((Ln(k) + LUff - 90)/10)).^alfaff - ..
                  (0.4 * 10.^((Tff + LUff - 90)/10)).^alfaff + ..
                  0.005135
            Lp(k,:) = 40*log10(Bff) + 94 
         end        
      else
         // Sound pressure level
         for k=1:length(Ln)
            Aff = 4.47e-3 * (10.^(0.025*Ln(k)) - 1.15) + ..
                  (0.4 * 10.^((Tff + LUff - 90)/10)).^alfaff;
            Lp(k,:) = 10./alfaff .* log10(Aff) - LUff + 94;
         end
      end
   
      
   case 2023
      f23  = [      20    25  31.5    40    50    63    80   100 ..
                   125   160   200   250   315   400   500   630 ..
                   800  1000  1250  1600  2000  2500  3150  4000 ..
                  5000  6300  8000 10000 12500 ];
      alfaf23 = [0.635 0.605 0.568 0.537 0.508 0.480 0.455 0.432 ..
                 0.410 0.389 0.371 0.356 0.340 0.329 0.318 0.310 ..
                 0.305 0.300 0.296 0.292 0.289 0.288 0.288 0.288 ..
                 0.289 0.294 0.303 0.320 0.354 ];
      LUf23 = [  -31.6 -27.2 -23.0 -19.1 -15.9 -13.0 -10.3  -8.1 ..
                  -6.2  -4.5  -3.1  -2.0  -1.1  -0.4   0.0   0.3 ..
                   0.5   0.0  -2.7  -4.1  -1.0   1.7   2.5   1.2 ..
                  -2.1  -7.1 -11.2 -10.7  -3.1]; 
      Tf23 = [    78.1  68.7  59.5  51.1  44.0  37.5  31.5  26.5 ..
                  22.1  17.9  14.4  11.4   8.6   6.2   4.4   3.0 ..
                   2.2   2.4   3.5   1.7  -1.3  -4.2  -6.0  -5.4 ..
                  -1.5   6.0  12.6  13.9  12.3];
      
      // Reference values
      alfar23 = 0.300;
      Tr23 = 2.4;
      po = 20e-6;
      
      // Interpolation of the parameters alfaf, 
      // LUf, Tf at f
      alfaff = interp1(f23, alfaf23, f, "spline");
      LUff = interp1(f23, LUf23, f, "spline");
      Tff = interp1(f23, Tf23, f, "spline");
      Lt = Tff;
      
      if di==-1
         // Loudness level
         for k=1:length(Ln)
            P = 10.^(alfaff  .* (Ln(k)+LUff)/10) - ..
                10.^(alfar23 .* (Tff + LUff)/10);
            Lp(k,:) = 10/alfar23 * log10(P.*(po^2).^(alfaff-alfar23) + ..
                      10.^(alfar23*Tr23/10));
         end
      else
         // Sound pressure level
         for k=1:length(Ln)
            Lp(k,:) = 10./alfaff .* log10((po^2).^(alfar23 - alfaff) .* ..
                  (10.^(alfar23*Ln(k)/10) - 10.^(alfar23*Tr23/10)) + ..
                  10.^(alfar23*(Tff + LUff)/10)) - LUff;
         end               
      end           
   end
   
endfunction


// Appendix
// The following values of cf61, which appear in
// the original Robinson & Dadson paper on which
// ISO R 266 is based, are seemingly incorrect. 
// They appear to have been corrected in the ISO
// recommendation, yielding smoother contours 
// for high loudness levels and low frequencies
//
//cf61 = [  -927   -782   -677   -589  -514  -454  -408  -371 ..
//          -343   -305   -285   -275  -268  -260  -255  -249 ..
//          -245   -241   -230   -218  -204  -188  -151  -116 ..
//           -84    -48    -22      0    23    47    71    90 ..  
//           100    118    120    113    98    80    59    14 ..
//           -35    -89   -145   -203  -269  -350  -454  -591 ..
//          -772]/1e5;
