text:
/*
 * SpanDSP - a series of DSP components for telephony
 *
 * make_line_models.c
 *
 * Written by Steve Underwood <steveu@coppice.org>
 *
 * Copyright (C) 2004 Steve Underwood
 *
 * All rights reserved.
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License version 2, as
 * published by the Free Software Foundation.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*! \page make_line_models_page Telephony line model construction
\section make_line_models_page_sec_1 What does it do?
???.

\section make_line_models_page_sec_2 How does it work?
???.
*/

#if defined(HAVE_CONFIG_H)
#include "config.h"
#endif

#include <stdlib.h>
#include <unistd.h>
#include <inttypes.h>
#include <string.h>
#include <stdio.h>
#if defined(HAVE_FFTW3_H)
#include <fftw3.h>
#else
#include <fftw.h>
#endif
#include <math.h>
#if defined(HAVE_STDBOOL_H)
#include <stdbool.h>
#else
#include "spandsp/stdbool.h"
#endif

#include "spandsp.h"

#if !defined(M_PI)
# define M_PI           3.14159265358979323846  /* pi */
#endif

#define LINE_MODEL_FILE_NAME    "line_models.c"

#define SAMPLE_RATE             8000

#define LINE_FILTER_SIZE        129
#define FFT_SIZE                1024

/* Tabulated medium range telephone line response
   (from p 537, Digital Communication, John G. Proakis */
/*
    amp 1.0 -> 2.15, freq = 3000 Hz -> 3.2, by 0.2 increments
    delay = 4 ms -> 2.2
 */

struct
{
    int frequency;
    float amp;
    float delay;
} proakis[] =
{
    {   0, 0.00,  4.80},
    { 200, 0.90,  3.50},
    { 400, 1.40,  2.20},
    { 600, 1.80,  0.90},
    { 800, 2.00,  0.50},
    {1000, 2.10,  0.25},
    {1200, 2.30,  0.10},
    {1400, 2.30,  0.05},
    {1600, 2.20,  0.00},
    {1800, 2.10,  0.00},
    {2000, 2.00,  0.00},
    {2200, 1.85,  0.05},
    {2400, 1.75,  0.10},
    {2600, 1.55,  0.20},
    {2800, 1.30,  0.40},
    {3000, 1.10,  0.50},
    {3200, 0.80,  0.90},
    {3400, 0.55,  1.20},
    {3600, 0.25,  2.20},
    {3800, 0.05,  3.20},
    {4000, 0.05,  4.20},
    {4200, 0.05,  5.20}
};

#define CPE_TO_CO_ATTENUATION       0       /* In dB */
#define CPE_TO_CO_DELAY             1       /* In us */
#define CPE_TO_CO_IMPEDANCE         2       /* In ohms */
#define CPE_TO_CO_PHASE             3       /* In degrees */
#define CO_TO_CPE_IMPEDANCE         4       /* In ohms */
#define CO_TO_CPE_PHASE             5       /* In degrees */
#define CO_TO_CPE_ATTENUATION       6       /* In dB */
#define CO_TO_CPE_DELAY             7       /* In us */

/* Terms used, for V.56bis:
    AD = attenuation distortion
    EDD = envelope delay distortion */

/* V.56bis EIA LL-1, non-loaded loop */

struct
{
    int freq;
    float ad[8];
} eia_ll1[] =
{
    { 200,    {0.0,    0.4,    767,     -1.4,    767,     -1.4,    0.0,    0.4}},
    { 300,    {0.0,    0.7,    766,     -2.0,    766,     -2.0,    0.0,    0.7}},
    { 400,    {0.0,    0.5,    763,     -2.8,    763,     -2.8,    0.0,    0.5}},
    { 500,    {0.0,    0.6,    765,     -3.4,    765,     -3.4,    0.0,    0.6}},
    { 600,    {0.0,    0.2,    764,     -4.1,    764,     -4.1,    0.0,    0.2}},
    { 700,    {0.0,    0.4,    764,     -4.7,    764,     -4.7,    0.0,    0.4}},
    { 800,    {0.0,    0.4,    762,     -5.4,    762,     -5.4,    0.0,    0.4}},
    { 900,    {0.0,    0.2,    762,     -6.0,    762,     -6.0,    0.0,    0.2}},
    {1000,    {1.2,    0.5,    761,     -6.7,    761,     -6.7,    1.2,    0.5}},
    {1100,    {0.0,    0.6,    759,     -7.4,    759,     -7.4,    0.0,    0.6}},
    {1200,    {0.0,    0.4,    757,     -8.1,    757,     -8.1,    0.0,    0.4}},
    {1300,    {0.0,    0.1,    757,     -8.6,    757,     -8.6,    0.0,    0.1}},
    {1400,    {0.0,    0.3,    755,     -9.3,    755,     -9.3,    0.0,    0.3}},
    {1500,    {0.0,    0.4,    753,    -10.0,    753,    -10.0,    0.0,    0.4}},
    {1600,    {0.0,    0.3,    751,    -10.7,    751,    -10.7,    0.0,    0.3}},
    {1700,    {0.0,    0.1,    748,    -11.3,    748,    -11.3,    0.0,    0.1}},
    {1800,    {0.0,   11.0,    748,    -11.9,    748,    -11.9,    0.0,   11.0}},
    {1900,    {0.1,    0.1,    745,    -12.5,    745,    -12.5,    0.1,    0.1}},
    {2000,    {0.1,    0.3,    743,    -13.9,    743,    -13.9,    0.1,    0.3}},
    {2100,    {0.1,    0.3,    740,    -13.9,    740,    -13.9,    0.1,    0.3}},
    {2200,    {0.1,    0.3,    737,    -14.5,    737,    -14.5,    0.1,    0.3}},
    {2300,    {0.1,    0.3,    734,    -15.2,    734,    -15.2,    0.1,    0.3}},
    {2400,    {0.1,    0.2,    731,    -15.8,    731,    -15.8,    0.1,    0.2}},
    {2500,    {0.1,    0.0,    728,    -16.4,    728,    -16.4,    0.1,    0.0}},
    {2600,    {0.1,    0.0,    729,    -16.8,    729,    -16.8,    0.1,    0.0}},
    {2700,    {0.2,    0.1,    726,    -17.4,    726,    -17.4,    0.2,    0.1}},
    {2800,    {0.2,    0.2,    722,    -18.0,    722,    -18.0,    0.2,    0.2}},
    {2900,    {0.2,    0.3,    719,    -18.6,    719,    -18.6,    0.2,    0.3}},
    {3000,    {0.2,    0.4,    715,    -19.3,    715,    -19.3,    0.2,    0.4}},
    {3100,    {0.2,    0.4,    712,    -19.9,    712,    -19.9,    0.2,    0.4}},
    {3200,    {0.2,    0.5,    708,    -20.5,    708,    -20.5,    0.2,    0.5}},
    {3300,    {0.2,    0.5,    704,    -21.1,    704,    -21.1,    0.2,    0.5}},
    {3400,    {0.2,    0.5,    700,    -21.7,    700,    -21.7,    0.2,    0.5}},
    {3500,    {0.2,    0.5,    696,    -22.3,    696,    -22.3,    0.2,    0.5}},
    {3600,    {0.2,    0.4,    692,    -22.9,    692,    -22.9,    0.2,    0.4}},
    {3700,    {0.2,    0.3,    688,    -23.5,    688,    -23.5,    0.2,    0.3}},
    {3800,    {0.2,    0.2,    684,    -24.1,    684,    -24.1,    0.2,    0.2}},
    {3900,    {0.2,    0.1,    680,    -24.7,    680,    -24.7,    0.2,    0.1}},
    {4000,    {0.2,   -0.1,    676,    -25.2,    676,    -25.2,    0.2,   -0.1}}
};

/* V.56bis EIA LL-2, non-loaded loop */

struct
{
    int freq;
    float ad[8];
} eia_ll2[] =
{
    { 200,   {-0.2,    6.6,   1086,     -4.9,   1085,     -5.6,   -0.2,    6.6}},
    { 300,   {-0.2,    7.7,   1079,     -7.3,   1077,     -8.3,   -0.2,    7.7}},
    { 400,   {-0.2,    6.7,   1062,     -9.9,   1058,    -11.2,   -0.2,    6.7}},
    { 500,   {-0.2,    7.1,   1059,    -12.0,   1053,    -13.6,   -0.2,    7.1}},
    { 600,   {-0.1,    5.2,   1041,    -14.4,   1034,    -16.3,   -0.1,    5.2}},
    { 700,   {-0.1,    5.8,   1030,    -16.5,   1020,    -18.6,   -0.1,    5.8}},
    { 800,   {-0.1,    5.4,   1010,    -18.7,    998,    -21.0,   -0.1,    5.4}},
    { 900,   { 0.0,    4.5,    997,    -20.5,    982,    -23.1,    0.0,    4.5}},
    {1000,   { 3.2,    5.1,    976,    -22.5,    959,    -25.3,    3.2,    5.1}},
    {1100,   { 0.1,    5.0,    954,    -24.5,    934,    -27.4,    0.1,    5.0}},
    {1200,   { 0.1,    4.0,    931,    -26.2,    909,    -29.4,    0.1,    4.0}},
    {1300,   { 0.2,    2.7,    918,    -27.6,    894,    -30.9,    0.2,    2.7}},
    {1400,   { 0.2,    2.8,    897,    -29.2,    871,    -32.6,    0.2,    2.8}},
    {1500,   { 0.3,    2.6,    874,    -30.7,    847,    -34.3,    0.3,    2.6}},
    {1600,   { 0.3,    2.0,    852,    -32.1,    823,    -35.8,    0.3,    2.0}},
    {1700,   { 0.4,    0.9,    831,    -33.4,    800,    -37.2,    0.4,    0.9}},
    {1800,   { 0.5,   40.8,    816,    -34.4,    783,    -38.4,    0.5,   40.8}},
    {1900,   { 0.6,    0.0,    796,    -35.6,    762,    -39.6,    0.6,    0.0}},
    {2000,   { 0.6,   -0.2,    776,    -36.6,    741,    -40.7,    0.6,   -0.2}},
    {2100,   { 0.7,   -0.6,    756,    -37.6,    720,    -41.8,    0.7,   -0.6}},
    {2200,   { 0.8,   -1.1,    737,    -38.6,    700,    -42.9,    0.8,   -1.1}},
    {2300,   { 0.9,   -1.8,    719,    -39.4,    681,    -43.8,    0.9,   -1.8}},
    {2400,   { 1.0,   -2.6,    701,    -40.2,    663,    -44.7,    1.0,   -2.6}},
    {2500,   { 1.0,   -3.7,    684,    -41.0,    646,    -45.5,    1.0,   -3.7}},
    {2600,   { 1.1,   -4.1,    678,    -41.3,    639,    -45.9,    1.1,   -4.1}},
    {2700,   { 1.2,   -4.3,    663,    -42.0,    623,    -46.6,    1.2,   -4.3}},
    {2800,   { 1.3,   -4.5,    647,    -42.6,    607,    -47.3,    1.3,   -4.5}},
    {2900,   { 1.4,   -4.8,    632,    -43.1,    591,    -47.9,    1.4,   -4.8}},
    {3000,   { 1.5,   -5.2,    617,    -43.6,    576,    -48.4,    1.5,   -5.2}},
    {3100,   { 1.6,   -5.6,    603,    -44.1,    562,    -49.0,    1.6,   -5.6}},
    {3200,   { 1.7,   -6.0,    589,    -44.5,    548,    -49.5,    1.7,   -6.0}},
    {3300,   { 1.8,   -6.5,    576,    -44.9,    535,    -49.9,    1.8,   -6.5}},
    {3400,   { 1.9,   -7.1,    563,    -45.3,    522,    -50.3,    1.9,   -7.1}},
    {3500,   { 2.0,   -7.7,    551,    -45.6,    509,    -50.7,    2.0,   -7.7}},
    {3600,   { 2.1,   -8.4,    539,    -45.9,    498,    -51.0,    2.1,   -8.4}},
    {3700,   { 2.2,   -9.1,    528,    -46.2,    486,    -51.3,    2.2,   -9.1}},
    {3800,   { 2.3,   -9.9,    518,    -46.4,    476,    -51.6,    2.3,   -9.9}},
    {3900,   { 2.4,  -10.6,    507,    -46.6,    466,    -51.9,    2.4,  -10.6}},
    {4000,   { 2.5,  -11.5,    498,    -46.8,    456,    -52.1,    2.5,  -11.5}}
};

/* V.56bis EIA LL-3, non-loaded loop */

struct
{
    int freq;
    float ad[8];
} eia_ll3[] =
{
    { 200,   {-0.3,   10.5,   1176,     -5.9,   1173,     -7.4,   -0.3,   10.5}},
    { 300,   {-0.3,   11.5,   1165,     -8.8,   1159,    -11.0,   -0.3,   11.5}},
    { 400,   {-0.3,   10.6,   1140,    -11.8,   1130,    -14.7,   -0.3,   10.6}},
    { 500,   {-0.3,   11.0,   1133,    -14.3,   1117,    -17.8,   -0.3,   11.0}},
    { 600,   {-0.2,    8.5,   1108,    -17.1,   1086,    -21.2,   -0.2,    8.5}},
    { 700,   {-0.2,    8.5,   1090,    -19.4,   1062,    -24.0,   -0.2,    8.5}},
    { 800,   {-0.1,    8.4,   1062,    -21.9,   1029,    -27.0,   -0.1,    8.4}},
    { 900,   { 0.0,    7.1,   1042,    -23.9,   1003,    -29.4,    0.0,    7.1}},
    {1000,   { 3.8,    7.7,   1013,    -23.0,    969,    -31.9,    3.8,    7.7}},
    {1100,   { 0.1,    7.4,    982,    -28.1,    934,    -34.3,    0.1,    7.4}},
    {1200,   { 0.1,    6.0,    953,    -29.9,    900,    -36.5,    0.1,    6.0}},
    {1300,   { 0.2,    4.2,    935,    -31.3,    878,    -38.1,    0.2,    4.2}},
    {1400,   { 0.3,    4.2,    907,    -32.8,    847,    -40.0,    0.3,    4.2}},
    {1500,   { 0.4,    3.7,    880,    -34.3,    817,    -41.7,    0.4,    3.7}},
    {1600,   { 0.5,    2.7,    853,    -35.6,    787,    -43.2,    0.5,    2.7}},
    {1700,   { 0.6,    1.2,    827,    -36.8,    760,    -44.6,    0.6,    1.2}},
    {1800,   { 0.7,   48.7,    809,    -37.8,    739,    -45.8,    0.7,   48.7}},
    {1900,   { 0.8,   -0.2,    785,    -38.8,    715,    -47.0,    0.8,   -0.2}},
    {2000,   { 0.9,   -0.7,    763,    -39.7,    691,    -48.0,    0.9,   -0.7}},
    {2100,   { 1.0,   -1.3,    741,    -40.5,    668,    -49.1,    1.0,   -1.3}},
    {2200,   { 1.1,   -2.1,    719,    -41.3,    647,    -50.0,    1.1,   -2.1}},
    {2300,   { 1.2,   -2.1,    699,    -42.0,    625,    -50.8,    1.2,   -2.1}},
    {2400,   { 1.2,   -4.3,    680,    -42.6,    606,    -51.6,    1.2,   -4.3}},
    {2500,   { 1.3,   -5.6,    663,    -43.2,    588,    -52.3,    1.3,   -5.6}},
    {2600,   { 1.6,   -6.2,    656,    -43.4,    581,    -52.7,    1.6,   -6.2}},
    {2700,   { 1.7,   -6.6,    640,    -43.9,    564,    -53.3,    1.7,   -6.6}},
    {2800,   { 1.8,   -7.0,    624,    -44.3,    548,    -53.9,    1.8,   -7.0}},
    {2900,   { 1.9,   -7.5,    609,    -44.7,    533,    -54.4,    1.9,   -7.5}},
    {3000,   { 2.0,   -8.0,    594,    -45.0,    518,    -54.8,    2.0,   -8.0}},
    {3100,   { 2.2,   -8.6,    580,    -45.3,    504,    -55.3,    2.2,   -8.6}},
    {3200,   { 2.3,   -9.2,    566,    -45.6,    490,    -55.7,    2.3,   -9.2}},
    {3300,   { 2.4,   -9.9,    553,    -45.8,    477,    -56.0,    2.4,   -9.9}},
    {3400,   { 2.6,  -10.7,    540,    -46.0,    465,    -56.3,    2.6,  -10.7}},
    {3500,   { 2.7,  -11.4,    529,    -46.2,    454,    -56.6,    2.7,  -11.4}},
    {3600,   { 2.8,  -12.3,    517,    -46.3,    443,    -56.9,    2.8,  -12.3}},
    {3700,   { 3.0,  -13.1,    507,    -46.4,    432,    -57.1,    3.0,  -13.1}},
    {3800,   { 3.1,  -14.0,    497,    -46.5,    422,    -57.3,    3.1,  -14.0}},
    {3900,   { 3.2,  -14.9,    487,    -46.6,    413,    -57.5,    3.2,  -14.9}},
    {4000,   { 3.3,  -15.9,    478,    -46.6,    404,    -57.7,    3.3,  -15.9}}
};

/* V.56bis EIA LL-4, non-loaded loop */

struct
{
    int freq;
    float ad[8];
} eia_ll4[] =
{
    { 200,   {-0.8,   31.0,   1564,    -10.7,   1564,    -10.7,   -0.8,   31.0}},
    { 300,   {-0.8,   32.6,   1520,    -15.6,   1520,    -15.6,   -0.8,   32.6}},
    { 400,   {-0.8,   29.8,   1447,    -20.5,   1447,    -20.5,   -0.8,   29.8}},
    { 500,   {-0.6,   29.7,   1402,    -24.3,   1402,    -24.3,   -0.6,   29.7}},
    { 600,   {-0.5,   24.9,   1328,    -28.1,   1328,    -28.1,   -0.5,   24.9}},
    { 700,   {-0.4,   24.8,   1270,    -31.2,   1270,    -31.2,   -0.4,   24.8}},
    { 800,   {-0.3,   22.7,   1200,    -34.0,   1200,    -34.0,   -0.3,   22.7}},
    { 900,   {-0.1,   19.8,   1148,    -36.2,   1148,    -36.2,   -0.1,   19.8}},
    {1000,   { 6.1,   19.3,   1086,    -38.3,   1086,    -38.3,    6.1,   19.3}},
    {1100,   { 0.1,   17.5,   1027,    -40.1,   1027,    -40.1,    0.1,   17.5}},
    {1200,   { 0.3,   14.3,    974,    -41.6,    974,    -41.6,    0.3,   14.3}},
    {1300,   { 0.5,   10.9,    941,    -42.6,    941,    -42.6,    0.5,   10.9}},
    {1400,   { 0.7,    9.6,    897,    -43.7,    897,    -43.7,    0.7,    9.6}},
    {1500,   { 0.9,    7.7,    856,    -44.6,    856,    -44.6,    0.9,    7.7}},
    {1600,   { 1.1,    5.3,    818,    -45.3,    818,    -45.3,    1.1,    5.3}},
    {1700,   { 1.3,    2.4,    784,    -45.9,    784,    -45.9,    1.3,    2.4}},
    {1800,   { 1.4,   69.1,    761,    -46.3,    761,    -46.3,    1.4,   69.1}},
    {1900,   { 1.7,   -1.3,    732,    -46.6,    732,    -46.6,    1.7,   -1.3}},
    {2000,   { 1.9,   -2.7,    706,    -46.9,    706,    -46.9,    1.9,   -2.7}},
    {2100,   { 2.1,   -4.3,    682,    -47.1,    682,    -47.1,    2.1,   -4.3}},
    {2200,   { 2.3,   -6.0,    659,    -47.3,    659,    -47.3,    2.3,   -6.0}},
    {2300,   { 2.5,   -7.9,    638,    -47.4,    638,    -47.4,    2.5,   -7.9}},
    {2400,   { 2.7,   -9.9,    619,    -47.4,    619,    -47.4,    2.7,   -9.9}},
    {2500,   { 2.9,  -12.0,    602,    -47.5,    602,    -47.5,    2.9,  -12.0}},
    {2600,   { 3.1,  -13.0,    596,    -47.4,    596,    -47.4,    3.1,  -13.0}},
    {2700,   { 3.3,  -13.9,    580,    -47.4,    580,    -47.4,    3.3,  -13.9}},
    {2800,   { 3.5,  -14.8,    566,    -47.3,    566,    -47.3,    3.5,  -14.8}},
    {2900,   { 3.7,  -15.7,    552,    -47.2,    552,    -47.2,    3.7,  -15.7}},
    {3000,   { 3.9,  -16.7,    539,    -47.1,    539,    -47.1,    3.9,  -16.7}},
    {3100,   { 4.1,  -17.7,    526,    -47.0,    526,    -47.0,    4.1,  -17.7}},
    {3200,   { 4.3,  -18.7,    515,    -46.8,    515,    -46.8,    4.3,  -18.7}},
    {3300,   { 4.5,  -19.8,    504,    -46.7,    504,    -46.7,    4.5,  -19.8}},
    {3400,   { 4.7,  -20.8,    493,    -46.5,    493,    -46.5,    4.7,  -20.8}},
    {3500,   { 4.9,  -21.8,    484,    -46.4,    484,    -46.4,    4.9,  -21.8}},
    {3600,   { 5.1,  -22.9,    475,    -46.2,    475,    -46.2,    5.1,  -22.9}},
    {3700,   { 5.3,  -23.9,    466,    -46.0,    466,    -46.0,    5.3,  -23.9}},
    {3800,   { 5.5,  -25.0,    458,    -45.9,    458,    -45.9,    5.5,  -25.0}},
    {3900,   { 5.6,  -26.1,    451,    -45.7,    451,    -45.7,    5.6,  -26.1}},
    {4000,   { 5.8,  -27.2,    444,    -45.5,    444,    -45.5,    5.8,  -27.2}}
};

/* V.56bis EIA LL-5, non-loaded loop */

struct
{
    int freq;
    float ad[8];
} eia_ll5[] =
{
    { 200,   {-1.4,    55.8,   1607,   -12.7,   1574,    -17.4,   -1.4,   55.8}},
    { 300,   {-1.3,    57.2,   1541,   -18.3,   1478,    -24.8,   -1.3,   57.2}},
    { 400,   {-1.2,    52.2,   1443,   -23.6,   1350,    -31.5,   -1.2,   52.2}},
    { 500,   {-1.0,    51.0,   1379,   -27.5,   1261,    -36.4,   -1.0,   51.0}},
    { 600,   {-0.9,    43.2,   1287,   -31.2,   1150,    -40.7,   -0.9,   43.2}},
    { 700,   {-0.7,    41.8,   1216,   -34.0,   1066,    -44.0,   -0.7,   41.8}},
    { 800,   {-0.5,    37.4,   1137,   -36.5,    979,    -46.9,   -0.5,   37.4}},
    { 900,   {-0.2,    32.4,   1080,   -38.3,    915,    -48.9,   -0.2,   32.4}},
    {1000,   { 7.0,    30.5,   1015,   -39.8,    848,    -50.7,    7.0,   30.5}},
    {1100,   { 0.3,    26.8,    956,   -41.1,    788,    -52.2,    0.3,   26.8}},
    {1200,   { 0.5,    21.5,    904,   -42.1,    736,    -53.3,    0.5,   21.5}},
    {1300,   { 0.8,    16.6,    873,   -42.7,    703,    -54.1,    0.8,   16.6}},
    {1400,   { 1.0,    14.1,    832,   -43.2,    663,    -54.8,    1.0,   14.1}},
    {1500,   { 1.3,    10.9,    795,   -43.7,    627,    -55.3,    1.3,   10.9}},
    {1600,   { 1.6,     7.3,    762,   -44.0,    595,    -55.7,    1.6,    7.3}},
    {1700,   { 1.9,     3.2,    733,   -44.2,    567,    -56.0,    1.9,    3.2}},
    {1800,   { 2.2,    81.5,    713,   -44.3,    547,    -56.2,    2.2,   81.5}},
    {1900,   { 2.4,    -1.9,    689,   -44.4,    524,    -56.4,    2.4,   -1.9}},
    {2000,   { 2.7,    -3.9,    667,   -44.4,    503,    -56.5,    2.7,   -3.9}},
    {2100,   { 3.0,    -6.1,    646,   -44.4,    485,    -56.5,    3.0,   -6.1}},
    {2200,   { 3.3,    -8.3,    628,   -44.4,    466,    -56.5,    3.3,   -8.3}},
    {2300,   { 3.6,   -10.7,    610,   -44.3,    450,    -56.5,    3.6,  -10.7}},
    {2400,   { 3.8,   -13.1,    595,   -44.2,    436,    -56.4,    3.8,  -13.1}},
    {2500,   { 4.1,   -15.5,    581,   -44.1,    422,    -56.3,    4.1,  -15.5}},
    {2600,   { 4.3,   -16.7,    577,   -44.0,    417,    -56.2,    4.3,  -16.7}},
    {2700,   { 4.6,   -17.7,    565,   -43.9,    406,    -56.1,    4.6,  -17.7}},
    {2800,   { 4.8,   -18.8,    553,   -43.8,    395,    -56.0,    4.8,  -18.8}},
    {2900,   { 5.1,   -19.9,    542,   -43.7,    395,    -55.9,    5.1,  -19.9}},
    {3000,   { 5.4,   -21.0,    531,   -43.6,    375,    -55.7,    5.4,  -21.0}},
    {3100,   { 5.6,   -22.1,    521,   -43.5,    366,    -55.6,    5.6,  -22.1}},
    {3200,   { 5.9,   -23.2,    511,   -43.4,    357,    -55.4,    5.9,  -23.2}},
    {3300,   { 6.1,   -24.3,    502,   -43.3,    349,    -55.3,    6.1,  -24.3}},
    {3400,   { 6.4,   -25.4,    494,   -43.2,    341,    -55.1,    6.4,  -25.4}},
    {3500,   { 6.6,   -26.5,    486,   -43.1,    334,    -55.0,    6.6,  -26.5}},
    {3600,   { 6.9,   -27.6,    478,   -43.0,    327,    -54.8,    6.9,  -27.6}},
    {3700,   { 7.1,   -28.7,    471,   -42.9,    321,    -54.7,    7.1,  -28.7}},
    {3800,   { 7.3,   -29.9,    464,   -42.8,    315,    -54.6,    7.3,  -29.9}},
    {3900,   { 7.5,   -31.0,    458,   -42.7,    310,    -54.4,    7.5,  -31.0}},
    {4000,   { 7.8,   -32.1,    452,   -42.7,    304,    -54.3,    7.8,  -32.1}}
};

/* V.56bis EIA LL-6, non-loaded loop */

struct
{
    int freq;
    float ad[8];
} eia_ll6[] =
{
    { 200,   {-0.2,  -39.3,   1756,    -12.0,   1748,    -19.8,   -0.2,  -39.3}},
    { 300,   {-0.2,  -31.7,   1642,    -15.9,   1689,    -26.9,   -0.2,  -31.7}},
    { 400,   {-0.2,  -37.5,   1506,    -18.4,   1427,    -33.4,   -0.2,  -37.5}},
    { 500,   {-0.1,  -34.7,   1442,    -19.5,   1301,    -37.7,   -0.1,  -34.7}},
    { 600,   {-0.1,  -46.0,   1363,    -20.1,   1153,    -40.7,   -0.1,  -46.0}},
    { 700,   { 0.0,  -40.8,   1320,    -20.7,   1045,    -42.2,    0.0,  -40.8}},
    { 800,   { 0.0,  -40.1,   1269,    -21.5,    943,    -42.3,    0.0,  -40.1}},
    { 900,   { 0.0,  -40.6,   1227,    -22.5,    878,    -41.3,    0.0,  -40.6}},
    {1000,   { 6.6,  -28.0,   1161,    -23.4,    825,    -39.3,    6.6,  -28.0}},
    {1100,   { 0.0,  -16.5,   1082,    -23.5,    797,    -36.8,    0.0,  -16.5}},
    {1200,   {-0.1,    0.3,   1000,    -22.2,    798,    -34.4,   -0.1,    0.3}},
    {1300,   { 0.0,   -2.3,    943,    -19.3,    826,    -33.2,    0.0,   -2.3}},
    {1400,   { 0.0,   13.5,    896,    -14.0,    870,    -33.8,    0.0,   13.5}},
    {1500,   { 0.1,   22.6,    890,     -7.2,    916,    -36.8,    0.1,   22.6}},
    {1600,   { 0.3,   30.3,    940,     -0.3,    938,    -42.0,    0.3,   30.3}},
    {1700,   { 0.5,   12.5,   1052,      4.6,    929,    -48.0,    0.5,   12.5}},
    {1800,   { 0.8,  458.6,   1212,      6.9,    880,    -52.8,    0.8,  458.6}},
    {1900,   { 1.1,   -5.1,   1410,      3.5,    814,    -56.5,    1.1,   -5.1}},
    {2000,   { 1.4,   -5.0,   1579,     -3.6,    747,    -58.5,    1.4,   -5.0}},
    {2100,   { 1.5,    6.1,   1618,    -13.2,    688,    -58.8,    1.5,    6.1}},
    {2200,   { 1.5,   33.5,   1491,    -21.5,    646,    -57.7,    1.5,   33.5}},
    {2300,   { 1.4,   80.5,   1275,    -24.9,    625,    -55.6,    1.4,   80.5}},
    {2400,   { 1.3,  142.3,   1078,    -20.8,    633,    -53.8,    1.3,  142.3}},
    {2500,   { 1.4,  196.5,    985,     -9.3,    664,    -54.5,    1.4,  196.5}},
    {2600,   { 1.6,  214.5,   1045,      2.4,    692,    -57.6,    1.6,  214.5}},
    {2700,   { 2.4,  196.8,   1326,     13.7,    684,    -63.5,    2.4,  196.8}},
    {2800,   { 3.4,  150.4,   1887,     14.7,    637,    -68.3,    3.4,  150.4}},
    {2900,   { 4.3,  125.3,   2608,      1.3,    501,    -70.7,    4.3,  125.3}},
    {3000,   { 4.9,  174.6,   2730,    -21.8,    533,    -70.6,    4.9,  174.6}},
    {3100,   { 4.9,  380.0,   2094,    -33.7,    506,    -68.5,    4.9,  380.0}},
    {3200,   { 5.2,  759.3,   1642,    -21.3,    522,    -67.0,    5.2,  759.3}},
    {3300,   { 8.0,  680.1,   2348,      0.5,    531,    -72.9,    8.0,  680.1}},
    {3400,   {13.1,  237.8,   4510,    -20.9,    482,    -77.3,   13.1,  237.8}},
    {3500,   {18.2,  -18.8,   4116,    -59.6,    439,    -78.0,   18.2,  -18.8}},
    {3600,   {22.7, -145.4,   3041,    -74.4,    487,    -77.7,   22.7, -145.4}},
    {3700,   {26.8, -214.5,   2427,    -80.1,    383,    -77.1,   26.8, -214.5}},
    {3800,   {30.4, -257.0,   2054,    -82.7,    364,    -76.4,   30.4, -257.0}},
    {3900,   {33.7, -285.6,   1803,    -84.2,    348,    -75.0,   33.7, -285.6}},
    {4000,   {36.8, -306.2,   1621,    -85.1,    334,    -75.7,   36.8, -306.2}}
};

/* V.56bis EIA LL-7, non-loaded loop */

struct
{
    int freq;
    float ad[8];
} eia_ll7[] =
{
    { 200,   { 0.4,   -81.3,    1848,  -10.5,   1737,    -15.6,    0.4,  -81.3}},
    { 300,   { 0.3,   -68.9,    1785,  -16.2,   1585,    -21.6,    0.3,  -68.9}},
    { 400,   { 0.2,   -68.1,    1646,  -22.0,   1388,    -25.8,    0.2,  -68.1}},
    { 500,   { 0.1,   -57.0,    1528,  -26.2,   1247,    -27.7,    0.1,  -57.0}},
    { 600,   { 0.0,   -59.8,    1349,  -28.9,   1087,    -27.3,    0.0,  -59.8}},
    { 700,   { 0.0,   -45.0,    1205,  -29.1,    975,    -24.8,    0.0,  -45.0}},
    { 800,   {-0.1,   -36.9,    1064,  -26.8,    885,    -19.7,   -0.1,  -36.9}},
    { 900,   {-0.1,   -37.1,     989,  -22.6,    846,    -13.5,   -0.1,  -37.1}},
    {1000,   { 5.9,   -29.2,     944,  -16.6,    847,     -6.1,    5.9,  -29.2}},
    {1100,   { 0.1,   -30.8,     951,  -10.5,    900,      0.3,    0.1,  -30.8}},
    {1200,   { 0.2,   -40.7,    1008,   -5.9,    999,      4.9,    0.2,  -40.7}},
    {1300,   { 0.4,   -53.3,    1897,   -4.0,   1122,      4.6,    0.4,  -53.3}},
    {1400,   { 0.5,   -52.7,    1197,   -4.8,   1253,      1.9,    0.5,  -52.7}},
    {1500,   { 0.6,   -48.3,    1269,   -8.4,   1339,     -3.8,    0.6,  -48.3}},
    {1600,   { 0.6,   -38.0,    1274,  -13.2,   1337,    -10.4,    0.6,  -38.0}},
    {1700,   { 0.5,   -21.6,    1208,  -16.9,   1250,    -15.2,    0.5,  -21.6}},
    {1800,   { 0.4,   539.7,    1119,  -17.8,   1143,    -16.6,    0.4,  539.7}},
    {1900,   { 0.3,    35.4,    1027,  -14.7,   1036,    -13.7,    0.3,   35.4}},
    {2000,   { 0.3,    64.1,     989,   -7.9,    998,     -6.9,    0.3,   64.1}},
    {2100,   { 0.4,    76.1,    1045,    0.1,   1040,      1.0,    0.4,   76.1}},
    {2200,   { 0.6,    69.8,    1210,    5.3,   1197,      6.9,    0.6,   69.8}},
    {2300,   { 1.0,    55.9,    1460,    4.6,   1430,      5.4,    1.0,   55.9}},
    {2400,   { 1.2,    51.3,    1692,   -2.8,   1640,     -1.7,    1.2,   51.3}},
    {2500,   { 1.3,    72.6,    1730,  -13.4,   1666,    -11.5,    1.3,   72.6}},
    {2600,   { 1.3,   117.1,    1613,  -49.6,   1556,    -16.9,    1.3,  117.1}},
    {2700,   { 1.1,   222.5,    1371,  -19.5,   1334,    -16.1,    1.1,  222.5}},
    {2800,   { 1.1,   332.3,    1258,   -8.9,   1243,     -5.1,    1.1,  332.3}},
    {2900,   { 1.7,   356.1,    1474,    4.8,   1480,      8.4,    1.7,  356.1}},
    {3000,   { 2.8,   299.9,    2128,    6.6,   2143,      9.8,    2.8,  299.9}},
    {3100,   { 3.9,   309.4,    2813,  -10.5,   2882,     -7.1,    3.9,  309.4}},
    {3200,   { 4.4,   576.4,    2490,  -27.7,   2487,    -22.2,    4.4,  576.4}},
    {3300,   { 5.6,  1030.6,    2237,  -17.4,   2385,     -9.0,    5.6, 1030.6}},
    {3400,   {10.7,   570.2,    3882,  -19.2,   4855,    -14.9,   10.7,  570.2}},
    {3500,   {17.3,    83.5,    4116,  -57.4,   4649,    -63.5,   17.3,   83.5}},
    {3600,   {23.2,  -130.6,    3057,  -74.0,   3175,    -78.6,   23.2, -130.6}},
    {3700,   {28.3,  -153.9,    2432,  -80.0,   2471,    -83.1,   28.3, -153.9}},
    {3800,   {32.8,  -292.4,    2055,  -82.8,   2072,    -85.1,   32.8, -292.4}},
    {3900,   {36.9,  -249.9,    1803,  -84.2,   1811,    -86.1,   36.9, -249.9}},
    {4000,   {40.7,  -356.2,    1621,  -85.1,   1625,    -86.7,   40.7, -356.2}}
};

/* V.56bis ETSI LL-1, non-loaded loop */

struct
{
    int freq;
    float ad[8];
} etsi_ll1[] =
{
    { 200,  {-0.78,   14.0, 1248.5,     -9.7, 1248.5,     -9.7,  -0.78,   14.0}},
    { 300,  {-0.74,   10.0, 1220.9,    -14.3, 1220.9,    -14.3,  -0.74,   10.0}},
    { 400,  {-0.68,    8.0, 1185.2,    -18.6, 1185.2,    -18.6,  -0.68,    8.0}},
    { 500,  {-0.60,    7.0, 1143.9,    -22.6, 1143.9,    -22.6,  -0.60,    7.0}},
    { 600,  {-0.51,    6.0, 1099.0,    -26.2, 1099.0,    -26.2,  -0.51,    6.0}},
    { 700,  {-0.40,    5.6, 1052.5,    -29.5, 1052.5,    -29.5,  -0.40,    5.6}},
    { 800,  {-0.28,    5.3, 1005.9,    -32.4, 1005.9,    -32.4,  -0.28,    5.3}},
    { 900,  {-0.14,    5.0,  960.3,    -35.0,  960.3,    -35.0,  -0.14,    5.0}},
    {1000,  { 4.7,     4.6,  916.4,    -37.3,  916.4,    -37.3,   4.7,     4.6}},
    {1100,  { 0.16,    4.3,  874.6,    -39.3,  874.6,    -39.3,   0.16,    4.3}},
    {1200,  { 0.33,    3.6,  835.3,    -41.1,  835.3,    -41.1,   0.33,    3.6}},
    {1300,  { 0.49,    2.6,  798.5,    -42.6,  798.5,    -42.6,   0.49,    2.6}},
    {1400,  { 0.67,    2.0,  764.2,    -43.9,  764.2,    -43.9,   0.67,    2.0}},
    {1500,  { 0.85,    1.0,  732.3,    -45.1,  732.3,    -45.1,   0.85,    1.0}},
    {1600,  { 1.04,    0.6,  702.7,    -46.1,  702.7,    -46.1,   1.04,    0.6}},
    {1700,  { 1.23,    0.3,  675.3,    -47.0,  675.3,    -47.0,   1.23,    0.3}},
    {1800,  { 1.43,   40.0,  649.8,    -47.7,  649.8,    -47.7,   1.43,   40.0}},
    {1900,  { 1.63,   -1.0,  626.2,    -48.4,  626.2,    -48.4,   1.63,   -1.0}},
    {2000,  { 1.83,   -2.0,  604.3,    -48.9,  604.3,    -48.9,   1.83,   -2.0}},
    {2100,  { 2.03,   -3.3,  584.0,    -49.4,  584.0,    -49.4,   2.03,   -3.3}},
    {2200,  { 2.23,   -3.6,  565.1,    -49.8,  565.1,    -49.8,   2.23,   -3.6}},
    {2300,  { 2.44,   -4.3,  547.5,    -50.1,  547.5,    -50.1,   2.44,   -4.3}},
    {2400,  { 2.64,   -5.0,  531.1,    -50.4,  531.1,    -50.4,   2.64,   -5.0}},
    {2500,  { 2.84,   -6.1,  515.9,    -50.6,  515.9,    -50.6,   2.84,   -6.1}},
    {2600,  { 3.05,   -6.6,  501.6,    -50.8,  501.6,    -50.8,   3.05,   -6.6}},
    {2700,  { 3.25,   -7.3,  488.2,    -51.0,  488.2,    -51.0,   3.25,   -7.3}},
    {2800,  { 3.45,   -7.6,  475.7,    -51.1,  475.7,    -51.1,   3.45,   -7.6}},
    {2900,  { 3.65,   -8.3,  464.0,    -51.1,  464.0,    -51.1,   3.65,   -8.3}},
    {3000,  { 3.85,   -8.6,  453.0,    -51.2,  453.0,    -51.2,   3.85,   -8.6}},
    {3100,  { 4.04,   -9.3,  442.6,    -51.2,  442.6,    -51.2,   4.04,   -9.3}},
    {3200,  { 4.24,  -10.3,  432.9,    -51.2,  432.9,    -51.2,   4.24,  -10.3}},
    {3300,  { 4.43,  -10.6,  423.7,    -51.2,  423.7,    -51.2,   4.43,  -10.6}},
    {3400,  { 4.62,  -11.3,  415.1,    -51.2,  415.1,    -51.2,   4.62,  -11.3}},
    {3500,  { 4.81,  -11.6,  406.9,    -51.1,  406.9,    -51.1,   4.81,  -11.6}},
    {3600,  { 5.00,  -12.3,  399.1,    -51.1,  399.1,    -51.1,   5.00,  -12.3}},
    {3700,  { 5.19,  -13.0,  391.8,    -51.0,  391.8,    -51.0,   5.19,  -13.0}},
    {3800,  { 5.37,  -13.4,  384.9,    -51.0,  384.9,    -51.0,   5.37,  -13.4}},
    {3900,  { 5.56,  -13.8,  378.3,    -50.9,  378.3,    -50.9,   5.56,  -13.8}},
    {4000,  { 5.74,  -14.4,  372.0,    -50.8,  372.0,    -50.8,   5.74,  -14.4}}
};

/* V.56bis ETSI LL-2, non-loaded loop */

struct
{
    int freq;
    float ad[8];
} etsi_ll2[] =
{
    { 200,  {-0.10,   15.0,  850.3,     -3.4,  850.3,     -3.4,  -0.10,   15.0}},
    { 300,  {-0.09,    8.0,  848.1,     -5.1,  848.1,     -5.1,  -0.09,    8.0}},
    { 400,  {-0.09,    7.0,  845.1,     -6.7,  845.1,     -6.7,  -0.09,    7.0}},
    { 500,  {-0.08,    5.0,  841.3,     -8.4,  841.3,     -8.4,  -0.08,    5.0}},
    { 600,  {-0.07,    4.6,  836.7,    -10.0,  836.7,    -10.0,  -0.07,    4.6}},
    { 700,  {-0.06,    4.3,  831.3,    -11.6,  831.3,    -11.6,  -0.06,    4.3}},
    { 800,  {-0.04,    3.8,  825.3,    -13.2,  825.3,    -13.2,  -0.04,    3.8}},
    { 900,  {-0.02,    3.4,  818.6,    -14.8,  818.6,    -14.8,  -0.02,    3.4}},
    {1000,  { 1.80,    3.0,  811.4,    -16.3,  811.4,    -16.3,   1.8,     3.0}},
    {1100,  { 0.02,    2.6,  803.6,    -17.8,  803.6,    -17.8,   0.02,    2.6}},
    {1200,  { 0.04,    2.3,  795.3,    -19.3,  795.3,    -19.3,   0.04,    2.3}},
    {1300,  { 0.06,    1.3,  786.6,    -20.7,  786.6,    -20.7,   0.06,    1.3}},
    {1400,  { 0.09,    0.9,  777.5,    -22.1,  777.5,    -22.1,   0.09,    0.9}},
    {1500,  { 0.12,    0.6,  768.1,    -23.5,  768.1,    -23.5,   0.12,    0.6}},
    {1600,  { 0.15,    0.3,  758.4,    -24.8,  758.4,    -24.8,   0.15,    0.3}},
    {1700,  { 0.18,    0.0,  748.4,    -26.1,  748.4,    -26.1,   0.18,    0.0}},
    {1800,  { 0.21,     15,  738.4,    -27.3,  738.4,    -27.3,   0.21,   15.0}},
    {1900,  { 0.24,   -1.0,  728.1,    -28.5,  728.1,    -28.5,   0.24,   -1.0}},
    {2000,  { 0.28,   -2.3,  717.8,    -29.7,  717.8,    -29.7,   0.28,   -2.3}},
    {2100,  { 0.32,   -2.6,  707.4,    -30.8,  707.4,    -30.8,   0.32,   -2.6}},
    {2200,  { 0.36,   -3.0,  697.0,    -31.9,  697.0,    -31.9,   0.36,   -3.0}},
    {2300,  { 0.40,   -3.3,  686.6,    -33.0,  686.6,    -33.0,   0.40,   -3.3}},
    {2400,  { 0.44,   -3.6,  676.2,    -34.0,  676.2,    -34.0,   0.44,   -3.6}},
    {2500,  { 0.48,   -4.5,  665.9,    -35.0,  665.9,    -35.0,   0.48,   -4.5}},
    {2600,  { 0.53,   -5.4,  655.6,    -35.9,  655.6,    -35.9,   0.53,   -5.4}},
    {2700,  { 0.57,   -6.3,  645.5,    -36.8,  645.5,    -36.8,   0.57,   -6.3}},
    {2800,  { 0.62,   -6.6,  635.5,    -37.7,  635.5,    -37.7,   0.62,   -6.6}},
    {2900,  { 0.67,   -6.9,  625.6,    -38.6,  625.6,    -38.6,   0.67,   -6.9}},
    {3000,  { 0.72,   -7.5,  615.8,    -39.4,  615.8,    -39.4,   0.72,   -7.5}},
    {3100,  { 0.77,   -8.3,  606.2,    -40.2,  606.2,    -40.2,   0.77,   -8.3}},
    {3200,  { 0.82,   -8.6,  596.7,    -40.9,  596.7,    -40.9,   0.82,   -8.6}},
    {3300,  { 0.87,   -9.3,  587.4,    -41.6,  587.4,    -41.6,   0.87,   -9.3}},
    {3400,  { 0.92,   -9.6,  578.3,    -42.3,  578.3,    -42.3,   0.92,   -9.6}},
    {3500,  { 0.98,  -10.3,  569.3,    -43.0,  569.3,    -43.0,   0.98,  -10.3}},
    {3600,  { 1.03,  -10.6,  560.6,    -43.7,  560.6,    -43.7,   1.03,  -10.6}},
    {3700,  { 1.09,  -11.3,  552.0,    -44.3,  552.0,    -44.3,   1.09,  -11.3}},
    {3800,  { 1.14,  -11.6,  543.5,    -44.9,  543.5,    -44.9,   1.14,  -11.6}},
    {3900,  { 1.20,  -12.3,  535.3,    -45.4,  535.3,    -45.4,   1.20,  -12.3}},
    {4000,  { 1.26,  -13.3,  527.2,    -46.0,  527.2,    -46.0,   1.26,  -13.3}}
};

/* V.56bis AD-1 AD-5 AD-6 AD-7 AD-8 AD-9 */

struct
{
    int freq;
    float ad[6];
} ad[] =
{
    {   0,  {90.0,   90.0,   90.0,   90.0,   90.0,   90.0}},
    { 200,  { 6.0,    3.2,    3.0,    2.9,   11.6,   23.3}},
    { 300,  { 1.3,    1.4,    1.2,    1.1,    6.9,   13.9}},
    { 400,  { 0.0,    0.4,    0.3,    0.3,    4.0,    7.9}},
    { 500,  { 0.0,   -0.1,    0.0,    0.1,    2.0,    4.1}},
    { 600,  { 0.0,   -0.1,    0.0,    0.1,    1.2,    2.4}},
    { 700,  { 0.0,    0.1,    0.0,    0.0,    0.8,    1.7}},
    { 800,  { 0.0,    0.0,    0.0,   -0.1,    0.5,    1.1}},
    { 900,  { 0.0,    0.0,    0.0,   -0.1,    0.2,    0.4}},
    {1000,  { 0.0,    0.0,    0.0,    0.0,    0.0,    0.0}},
    {1100,  { 0.0,    0.0,    0.1,    0.0,   -0.1,   -0.2}},
    {1200,  { 0.0,    0.0,    0.1,    0.1,   -0.1,   -0.2}},
    {1300,  { 0.0,    0.1,    0.2,    0.3,   -0.1,   -0.2}},
    {1400,  { 0.0,    0.2,    0.3,    0.4,   -0.1,   -0.3}},
    {1500,  { 0.0,    0.2,    0.3,    0.4,   -0.2,   -0.4}},
    {1600,  { 0.0,    0.3,    0.5,    0.5,   -0.1,   -0.3}},
    {1700,  { 0.0,    0.3,    0.5,    0.6,   -0.1,   -0.1}},
    {1800,  { 0.0,    0.3,    0.5,    0.6,    0.0,    0.0}},
    {1900,  { 0.0,    0.4,    0.7,    0.7,    0.1,    0.2}},
    {2000,  { 0.0,    0.5,    0.8,    0.9,    0.2,    0.5}},
    {2100,  { 0.1,    0.6,    1.0,    1.0,    0.5,    0.9}},
    {2200,  { 0.2,    0.7,    1.1,    1.1,    0.6,    1.1}},
    {2300,  { 0.3,    0.9,    1.2,    1.4,    0.8,    1.5}},
    {2400,  { 0.4,    1.1,    1.5,    1.6,    0.9,    1.8}},
    {2500,  { 0.5,    1.3,    1.8,    2.0,    1.1,    2.3}},
    {2600,  { 0.6,    1.6,    2.4,    2.7,    1.4,    2.8}},
    {2700,  { 0.7,    2.0,    3.0,    3.5,    1.7,    3.4}},
    {2800,  { 0.7,    2.3,    3.5,    4.3,    2.0,    4.0}},
    {2900,  { 0.9,    2.8,    4.2,    5.0,    2.4,    4.9}},
    {3000,  { 1.1,    3.2,    4.9,    5.8,    3.0,    5.9}},
    {3100,  { 1.2,    3.5,    5.6,    6.7,    3.4,    6.8}},
    {3200,  { 1.3,    4.1,    6.7,    8.0,    3.9,    7.7}},
    {3300,  { 1.6,    4.8,    8.0,    9.6,    4.6,    9.2}},
    {3400,  { 1.8,    5.3,    9.1,   11.0,    5.4,   10.7}},
    {3500,  { 2.4,    5.7,   10.3,   12.2,    6.3,   12.6}},
    {3600,  { 3.0,    6.6,   12.1,   13.9,    7.8,   15.5}},
    {3700,  { 5.7,    8.9,   15.8,   17.3,   10.3,   20.5}},
    {3800,  {13.5,   15.7,   24.4,   25.7,   16.2,   32.4}},
    {3900,  {31.2,   31.1,   42.2,   43.3,   29.9,   59.9}},
    {4000,  {31.2,   31.1,   42.2,   43.3,   29.9,   59.9}}
};

/* V.56bis EDD-1 EDD-2 EDD-3 */

struct
{
    int freq;
    float edd[3];
} edd[] =
{
    {   0,    {3.98,    3.76,    8.00}},
    { 200,    {3.98,    3.76,    8.00}},
    { 300,    {2.70,    3.76,    8.00}},
    { 400,    {1.69,    2.20,    6.90}},
    { 500,    {1.15,    1.36,    5.50}},
    { 600,    {0.80,    0.91,    4.40}},
    { 700,    {0.60,    0.64,    3.40}},
    { 800,    {0.50,    0.46,    2.80}},
    { 900,    {0.40,    0.34,    2.00}},
    {1000,    {0.30,    0.24,    1.50}},
    {1100,    {0.20,    0.16,    1.00}},
    {1200,    {0.20,    0.11,    0.70}},
    {1300,    {0.10,    0.07,    0.40}},
    {1400,    {0.05,    0.05,    0.30}},
    {1500,    {0.00,    0.03,    0.20}},
    {1600,    {0.00,    0.01,    0.10}},
    {1700,    {0.00,    0.00,    0.10}},
    {1800,    {0.00,    0.00,    0.00}},
    {1900,    {0.00,    0.02,    0.10}},
    {2000,    {0.00,    0.04,    0.10}},
    {2100,    {0.02,    0.08,    0.10}},
    {2200,    {0.02,    0.12,    0.20}},
    {2300,    {0.02,    0.16,    0.20}},
    {2400,    {0.02,    0.20,    0.30}},
    {2500,    {0.10,    0.27,    0.40}},
    {2600,    {0.12,    0.36,    0.50}},
    {2700,    {0.15,    0.47,    0.80}},
    {2800,    {0.20,    0.60,    1.10}},
    {2900,    {0.27,    0.77,    1.50}},
    {3000,    {0.40,    1.01,    2.00}},
    {3100,    {0.56,    1.32,    2.60}},
    {3200,    {0.83,    1.78,    3.20}},
    {3300,    {1.07,    1.78,    4.00}},
    {3400,    {1.39,    1.78,    4.00}},
    {3500,    {1.39,    1.78,    4.00}},
    {3600,    {1.39,    1.78,    4.00}},
    {3700,    {1.39,    1.78,    4.00}},
    {3800,    {1.39,    1.78,    4.00}},
    {3900,    {1.39,    1.78,    4.00}},
    {4000,    {1.39,    1.78,    4.00}}
};

/* V.56bis PCM AD-1, AD-2, AD-3 */

struct
{
    int freq;
    float ad[3];
} pcm_ad[] =
{
    {  50,    {41.4,    77.8,   114.2}},
    { 100,    {15.5,    27.7,    39.9}},
    { 150,    { 3.7,     6.1,     8.6}},
    { 200,    { 0.5,     0.8,     1.0}},
    { 250,    {-0.2,    -0.2,    -0.3}},
    { 300,    {-0.2,    -0.3,    -0.4}},
    { 400,    { 0.0,    -0.2,    -0.3}},
    { 500,    {-0.2,    -0.4,    -0.5}},
    { 600,    {-0.2,    -0.3,    -0.5}},
    { 700,    {-0.2,    -0.3,    -0.5}},
    { 800,    {-0.2,    -0.4,    -0.5}},
    { 900,    {-0.2,    -0.3,    -0.4}},
    {1000,    {-0.1,    -0.2,    -0.3}},
    {1100,    {-0.2,    -0.3,    -0.3}},
    {1200,    {-0.2,    -0.3,    -0.4}},
    {1300,    {-0.2,    -0.3,    -0.5}},
    {1400,    {-0.1,    -0.3,    -0.4}},
    {1500,    {-0.1,    -0.3,    -0.4}},
    {1600,    {-0.1,    -0.2,    -0.3}},
    {1700,    {-0.1,    -0.3,    -0.4}},
    {1800,    {-0.2,    -0.3,    -0.4}},
    {1900,    {-0.2,    -0.3,    -0.3}},
    {2000,    {-0.1,    -0.2,    -0.3}},
    {2100,    {-0.1,    -0.2,    -0.3}},
    {2200,    {-0.1,    -0.3,    -0.4}},
    {2300,    {-0.1,    -0.1,    -0.2}},
    {2400,    {-0.1,    -0.1,    -0.2}},
    {2500,    { 0.0,    -0.1,    -0.1}},
    {2600,    { 0.0,    -0.1,    -0.1}},
    {2700,    { 0.0,     0.0,     0.1}},
    {2800,    { 0.0,     0.0,     0.1}},
    {2900,    { 0.1,     0.2,     0.2}},
    {3000,    { 0.0,     0.0,     0.1}},
    {3100,    { 0.0,     0.0,     0.0}},
    {3200,    { 0.0,     0.0,     0.1}},
    {3300,    { 0.3,     0.7,     1.0}},
    {3400,    { 1.2,     2.4,     3.6}},
    {3500,    { 3.2,     6.3,     9.5}},
    {3550,    { 5.0,     9.6,    14.3}},
    {3600,    { 7.0,    13.5,    19.9}},
    {3650,    {10.0,    18.7,    27.5}},
    {3700,    {13.4,    24.6,    35.8}},
    {3750,    {18.1,    32.1,    46.2}},
    {3800,    {24.3,    41.2,    58.2}},
    {3850,    {32.5,    52.6,    72.7}},
    {3900,    {43.4,    66.6,    89.8}},
    {4000,    {43.4,    66.6,    89.8}}
};

/* V.56bis PCM EDD-1, EDD-2, EDD-3 */

struct
{
    int freq;
    float edd[3];
} pcm_edd[] =
{
    { 150,   { 2.76,    5.5,    8.3}},
    { 200,   { 1.70,    3.4,    5.1}},
    { 250,   { 0.92,    1.8,    2.8}},
    { 300,   { 0.55,    1.1,    1.7}},
    { 400,   { 0.25,    0.5,    0.7}},
    { 500,   { 0.12,    0.2,    0.4}},
    { 600,   { 0.06,    0.1,    0.2}},
    { 700,   { 0.03,    0.1,    0.1}},
    { 800,   { 0.01,    0.0,    0.0}},
    { 900,   { 0.00,    0.0,    0.0}},
    {1000,   {-0.01,    0.0,    0.0}},
    {1100,   {-0.01,    0.0,    0.0}},
    {1200,   {-0.02,    0.0,   -0.1}},
    {1300,   {-0.02,    0.0,   -0.1}},
    {1400,   {-0.01,    0.0,    0.0}},
    {1500,   {-0.01,    0.0,    0.0}},
    {1600,   { 0.00,    0.0,    0.0}},
    {1700,   { 0.00,    0.0,    0.0}},
    {1800,   { 0.01,    0.0,    0.0}},
    {1900,   { 0.02,    0.0,    0.0}},
    {2000,   { 0.02,    0.0,    0.1}},
    {2100,   { 0.04,    0.1,    0.1}},
    {2200,   { 0.05,    0.1,    0.2}},
    {2300,   { 0.06,    0.1,    0.2}},
    {2400,   { 0.07,    0.1,    0.2}},
    {2500,   { 0.10,    0.2,    0.3}},
    {2600,   { 0.11,    0.2,    0.3}},
    {2700,   { 0.14,    0.3,    0.4}},
    {2800,   { 0.18,    0.4,    0.5}},
    {2900,   { 0.22,    0.4,    0.6}},
    {3000,   { 0.27,    0.5,    0.8}},
    {3100,   { 0.34,    0.7,    1.0}},
    {3200,   { 0.45,    0.9,    1.4}},
    {3250,   { 0.52,    1.0,    1.6}},
    {3300,   { 0.60,    1.2,    1.8}},
    {3350,   { 0.66,    1.3,    2.0}},
    {3400,   { 0.74,    1.5,    2.2}},
    {3450,   { 0.79,    1.6,    2.4}},
    {3500,   { 0.83,    1.7,    2.5}},
    {3550,   { 0.84,    1.7,    2.5}},
    {3600,   { 0.81,    1.6,    2.4}},
    {3700,   { 0.81,    1.6,    2.4}},
    {3800,   { 0.81,    1.6,    2.4}},
    {3900,   { 0.81,    1.6,    2.4}},
    {4000,   { 0.81,    1.6,    2.4}}
};

FILE *outfile;

float impulse_responses[100][LINE_FILTER_SIZE];
int filter_sets = 0;

static void generate_ad_edd(void)
{
    float f;
    float offset;
    float amp;
    float phase;
    //float delay;
    float pw;
#if defined(HAVE_FFTW3_H)
    double in[FFT_SIZE][2];
    double out[FFT_SIZE][2];
#else
    fftw_complex in[FFT_SIZE];
    fftw_complex out[FFT_SIZE];
#endif
    fftw_plan p;
    int i;
    int j;
    int k;
    int l;

#if defined(HAVE_FFTW3_H)
    p = fftw_plan_dft_1d(FFT_SIZE, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
#else
    p = fftw_create_plan(FFT_SIZE, FFTW_BACKWARD, FFTW_ESTIMATE);
#endif
    for (j = 0;  j < 6;  j++)
    {
        for (k = 0;  k < 3;  k++)
        {
            for (i = 0;  i < FFT_SIZE;  i++)
            {
#if defined(HAVE_FFTW3_H)
                in[i][0] =
                in[i][1] = 0.0f;
#else
                in[i].re =
                in[i].im = 0.0f;
#endif
            }
            for (i = 1;  i < FFT_SIZE/2;  i++)
            {
                f = (float) i*SAMPLE_RATE/FFT_SIZE;
                amp = 0.0f;
                for (l = 0;  l < (int) (sizeof(ad)/sizeof(ad[0]));  l++)
                {
                    if (f < ad[l].freq)
                        break;
                }
                if (l < (int) (sizeof(ad)/sizeof(ad[0])))
                {
                    offset = (f - ad[l - 1].freq)/(ad[l].freq - ad[l - 1].freq);
                    amp = (1.0 - offset)*ad[l - 1].ad[j] + offset*ad[l].ad[j];
                    amp = pow(10.0, -amp/20.0);
                }
                //delay = 0.0f;
                for (l = 0;  l < (int) (sizeof(edd)/sizeof(edd[0]));  l++)
                {
                    if (f < edd[l].freq)
                        break;
                }
                if (l < (int) (sizeof(edd)/sizeof(edd[0])))
                {
                    offset = (f - edd[l - 1].freq)/(edd[l].freq - edd[l - 1].freq);
                    //delay = (1.0f - offset)*edd[l - 1].edd[k] + offset*edd[l].edd[k];
                }
                //phase = 2.0f*M_PI*f*delay*0.001f;
                phase = 0.0f;
#if defined(HAVE_FFTW3_H)
                in[i][0] = amp*cosf(phase);
                in[i][1] = amp*sinf(phase);
                in[FFT_SIZE - i][0] = in[i][0];
                in[FFT_SIZE - i][1] = -in[i][1];
#else
                in[i].re = amp*cosf(phase);
                in[i].im = amp*sinf(phase);
                in[FFT_SIZE - i].re = in[i].re;
                in[FFT_SIZE - i].im = -in[i].im;
#endif
            }
#if 0
            for (i = 0;  i < FFT_SIZE;  i++)
                fprintf(outfile, "%5d %15.5f,%15.5f\n", i, in[i].re, in[i].im);
#endif
#if defined(HAVE_FFTW3_H)
            fftw_execute(p);
#else
            fftw_one(p, in, out);
#endif

            fprintf(outfile, "/* V.56bis AD-%d, EDD%d */\n", (j == 0)  ?  1  :  j + 4, k + 1);

            fprintf(outfile, "const float ad_%d_edd_%d_model[] =\n", (j == 0)  ?  1  :  j + 4, k + 1);
            fprintf(outfile, "{\n");
            /* Normalise the filter's gain */
            pw = 0.0f;
            l = FFT_SIZE - (LINE_FILTER_SIZE - 1)/2;
            for (i = 0;  i < LINE_FILTER_SIZE;  i++)
            {
#if defined(HAVE_FFTW3_H)
                pw += out[l][0]*out[l][0];
#else
                pw += out[l].re*out[l].re;
#endif
                if (++l == FFT_SIZE)
                    l = 0;
            }
            pw = sqrt(pw);
            l = FFT_SIZE - (LINE_FILTER_SIZE - 1)/2;
            for (i = 0;  i < LINE_FILTER_SIZE;  i++)
            {
#if defined(HAVE_FFTW3_H)
                impulse_responses[filter_sets][i] = out[l][0]/pw;
#else
                impulse_responses[filter_sets][i] = out[l].re/pw;
#endif
                fprintf(outfile, "%15.5f,\n", impulse_responses[filter_sets][i]);
                if (++l == FFT_SIZE)
                    l = 0;
            }
            fprintf(outfile, "};\n\n");
            filter_sets++;
        }
    }
}

static void generate_proakis(void)
{
    float f;
    float f1;
    float offset;
    float amp;
    float phase;
    //float delay;
    float pw;
    int index;
    int i;
    int l;
#if defined(HAVE_FFTW3_H)
    double in[FFT_SIZE][2];
    double out[FFT_SIZE][2];
#else
    fftw_complex in[FFT_SIZE];
    fftw_complex out[FFT_SIZE];
#endif
    fftw_plan p;

#if defined(HAVE_FFTW3_H)
    p = fftw_plan_dft_1d(FFT_SIZE, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
#else
    p = fftw_create_plan(FFT_SIZE, FFTW_BACKWARD, FFTW_ESTIMATE);
#endif
    for (i = 0;  i < FFT_SIZE;  i++)
    {
#if defined(HAVE_FFTW3_H)
        in[i][0] =
        in[i][1] = 0.0f;
#else
        in[i].re =
        in[i].im = 0.0f;
#endif
    }
    for (i = 1;  i < FFT_SIZE/2;  i++)
    {
        f = (float) i*SAMPLE_RATE/FFT_SIZE;
        f1 = f/200.0f;
        offset = f1 - floor(f1);
        index = (int) floor(f1);

        /* Linear interpolation */
        amp = ((1.0f - offset)*proakis[index].amp + offset*proakis[index + 1].amp)/2.3f;
        //delay = (1.0f - offset)*proakis[index].delay + offset*proakis[index + 1].delay;
        //phase = 2.0f*M_PI*f*delay*0.001f;
        phase = 0.0f;
#if defined(HAVE_FFTW3_H)
        in[i][0] = amp*cosf(phase);
        in[i][1] = amp*sinf(phase);
        in[FFT_SIZE - i][0] = in[i][0];
        in[FFT_SIZE - i][1] = -in[i][1];
#else
        in[i].re = amp*cosf(phase);
        in[i].im = amp*sinf(phase);
        in[FFT_SIZE - i].re = in[i].re;
        in[FFT_SIZE - i].im = -in[i].im;
#endif
    }

#if defined(HAVE_FFTW3_H)
    fftw_execute(p);
#else
    fftw_one(p, in, out);
#endif

    fprintf(outfile, "/* Medium range telephone line response\n");
    fprintf(outfile, "   (from p 537, Digital Communication, John G. Proakis */\n");

    fprintf(outfile, "const float proakis_line_model[] =\n");
    fprintf(outfile, "{\n");
    /* Normalise the filter's gain */
    pw = 0.0f;
    l = FFT_SIZE - (LINE_FILTER_SIZE - 1)/2;
    for (i = 0;  i < LINE_FILTER_SIZE;  i++)
    {
#if defined(HAVE_FFTW3_H)
        pw += out[l][0]*out[l][0];
#else
        pw += out[l].re*out[l].re;
#endif
        if (++l == FFT_SIZE)
            l = 0;
    }
    pw = sqrt(pw);
    l = FFT_SIZE - (LINE_FILTER_SIZE - 1)/2;
    for (i = 0;  i < LINE_FILTER_SIZE;  i++)
    {
#if defined(HAVE_FFTW3_H)
        impulse_responses[filter_sets][i] = out[l][0]/pw;
#else
        impulse_responses[filter_sets][i] = out[l].re/pw;
#endif
        fprintf(outfile, "%15.5f,\n", impulse_responses[filter_sets][i]);
        if (++l == FFT_SIZE)
            l = 0;
    }
    fprintf(outfile, "};\n\n");
    filter_sets++;
}

int main(int argc, char *argv[])
{
    int i;
    int j;

    if ((outfile = fopen(LINE_MODEL_FILE_NAME, "w")) == NULL)
    {
        fprintf(stderr, "Failed to open %s\n", "line_model.txt");
        exit(2);
    }

    generate_proakis();
    generate_ad_edd();

    fclose(outfile);

    if (argc > 1)
    {
        for (i = 0;  i < LINE_FILTER_SIZE;  i++)
        {
            printf("%d, ", i);
            for (j = 0;  j < filter_sets;  j++)
            {
                printf("%15.5f, ", impulse_responses[j][i]);
            }
            printf("\n");
        }
    }
    return 0;
}

max_stars_repo_name: zengfanmao/mpds

max_issues_repo_name: zengfanmao/mpds

max_forks_repo_name: zengfanmao/mpds


text:
// ---------------------------------------------------------------------------
//
//  Author  : github.com/luncliff (luncliff@gmail.com)
//  License : CC BY 4.0
//
//  Note
//      Async I/O operation support for socket
//
// ---------------------------------------------------------------------------
#pragma once
// clang-format off
#ifdef USE_STATIC_LINK_MACRO // ignore macro declaration in static build
#   define _INTERFACE_
#   define _HIDDEN_
#else
#   if defined(_MSC_VER) // MSVC
#       define _HIDDEN_
#       ifdef _WINDLL
#           define _INTERFACE_ __declspec(dllexport)
#       else
#           define _INTERFACE_ __declspec(dllimport)
#       endif
#   elif defined(__GNUC__) || defined(__clang__)
#       define _INTERFACE_ __attribute__((visibility("default")))
#       define _HIDDEN_ __attribute__((visibility("hidden")))
#   else
#       error "unexpected compiler"
#   endif // compiler check
#endif
// clang-format on

#ifndef COROUTINE_NET_IO_H
#define COROUTINE_NET_IO_H
#include <coroutine/yield.hpp>

#include <chrono>
#include <gsl/gsl>

#if defined(_MSC_VER)
#include <WS2tcpip.h>
#include <WinSock2.h>
#include <ws2def.h>
#pragma comment(lib, "Ws2_32.lib")

using io_control_block = OVERLAPPED;

static constexpr bool is_winsock = true;
static constexpr bool is_netinet = false;
#else
#include <fcntl.h>
#include <netdb.h>
#include <netinet/in.h>
#include <sys/socket.h>
#include <unistd.h>

// - Note
// follow the definition of Windows `OVERLAPPED`
struct io_control_block
{
    uint64_t internal; // uint32_t errc; int32_t flag;
    uint64_t internal_high;
    union {
        struct
        {
            int32_t offset; // socklen_t addrlen;
            int32_t offset_high;
        };
        void* ptr;
    };
    int64_t handle; // int64_t sd;
};

static constexpr bool is_winsock = false;
static constexpr bool is_netinet = true;
#endif

// - Note
//      Even though this class violates C++ Core Guidelines,
//      it will be used for simplicity
union endpoint_t final {
    sockaddr_storage storage{};
    sockaddr addr;
    sockaddr_in in4;
    sockaddr_in6 in6;
};

using io_task_t = std::experimental::coroutine_handle<void>;
using buffer_view_t = gsl::span<gsl::byte>;

struct io_work_t : public io_control_block
{
    io_task_t task{};
    buffer_view_t buffer{};
    endpoint_t* ep{};

  public:
    _INTERFACE_ bool ready() const noexcept;
    _INTERFACE_ uint32_t error() const noexcept;
};
static_assert(sizeof(buffer_view_t) <= sizeof(void*) * 2);
static_assert(sizeof(io_work_t) <= 64);

class io_send_to final : public io_work_t
{
  public:
    _INTERFACE_ void suspend(io_task_t rh) noexcept(false);
    _INTERFACE_ int64_t resume() noexcept;

  public:
    auto await_ready() const noexcept
    {
        return this->ready();
    }
    void await_suspend(io_task_t rh) noexcept(false)
    {
        return this->suspend(rh);
    }
    auto await_resume() noexcept
    {
        return this->resume();
    }
};
static_assert(sizeof(io_send_to) == sizeof(io_work_t));

[[nodiscard]] _INTERFACE_ auto
    send_to(uint64_t sd, const sockaddr_in& remote, //
            buffer_view_t buffer, io_work_t& work) noexcept(false)
        -> io_send_to&;

[[nodiscard]] _INTERFACE_ //
    auto
    send_to(uint64_t sd, const sockaddr_in6& remote, //
            buffer_view_t buffer, io_work_t& work) noexcept(false)
        -> io_send_to&;

class io_recv_from final : public io_work_t
{
  public:
    _INTERFACE_ void suspend(io_task_t rh) noexcept(false);
    _INTERFACE_ int64_t resume() noexcept;

  public:
    auto await_ready() const noexcept
    {
        return this->ready();
    }
    void await_suspend(io_task_t rh) noexcept(false)
    {
        return this->suspend(rh);
    }
    auto await_resume() noexcept
    {
        return this->resume();
    }
};
static_assert(sizeof(io_recv_from) == sizeof(io_work_t));

[[nodiscard]] _INTERFACE_ //
    auto
    recv_from(uint64_t sd, sockaddr_in6& remote, //
              buffer_view_t buffer, io_work_t& work) noexcept(false)
        -> io_recv_from&;

[[nodiscard]] _INTERFACE_ //
    auto
    recv_from(uint64_t sd, sockaddr_in& remote, //
              buffer_view_t buffer, io_work_t& work) noexcept(false)
        -> io_recv_from&;

class io_send final : public io_work_t
{
  public:
    _INTERFACE_ void suspend(io_task_t rh) noexcept(false);
    _INTERFACE_ int64_t resume() noexcept;

  public:
    auto await_ready() const noexcept
    {
        return this->ready();
    }
    void await_suspend(io_task_t rh) noexcept(false)
    {
        return this->suspend(rh);
    }
    auto await_resume() noexcept
    {
        return this->resume();
    }
};
static_assert(sizeof(io_send) == sizeof(io_work_t));

[[nodiscard]] _INTERFACE_ //
    auto
    send_stream(uint64_t sd, buffer_view_t buffer, uint32_t flag,
                io_work_t& work) noexcept(false) -> io_send&;

class io_recv final : public io_work_t
{
  public:
    _INTERFACE_ void suspend(io_task_t rh) noexcept(false);
    _INTERFACE_ int64_t resume() noexcept;

  public:
    auto await_ready() const noexcept
    {
        return this->ready();
    }
    void await_suspend(io_task_t rh) noexcept(false)
    {
        return this->suspend(rh);
    }
    auto await_resume() noexcept
    {
        return this->resume();
    }
};
static_assert(sizeof(io_recv) == sizeof(io_work_t));

[[nodiscard]] _INTERFACE_ //
    auto
    recv_stream(uint64_t sd, buffer_view_t buffer, uint32_t flag,
                io_work_t& work) noexcept(false) -> io_recv&;

// - Note
//      This function is only non-windows platform.
//      Over windows api, it always yields nothing.
//
//      User must continue looping without break so that there is no leak of
//      event. Also, the library doesn't guarantee all coroutines(i/o tasks)
//      will be fetched at once. Therefore it is strongly recommended for user
//      to have a method to detect that coroutines are returned.
_INTERFACE_
auto wait_io_tasks(std::chrono::nanoseconds timeout) noexcept(false)
    -> coro::enumerable<io_task_t>;

_INTERFACE_
auto host_name() noexcept -> gsl::czstring<NI_MAXHOST>;

_INTERFACE_
std::errc peer_name(uint64_t sd, sockaddr_in6& ep) noexcept;

_INTERFACE_
std::errc sock_name(uint64_t sd, sockaddr_in6& ep) noexcept;

_INTERFACE_
auto resolve(const addrinfo& hint, //
             gsl::czstring<NI_MAXHOST> name,
             gsl::czstring<NI_MAXSERV> serv) noexcept
    -> coro::enumerable<sockaddr_in6>;

_INTERFACE_
std::errc nameof(const sockaddr_in& ep, //
                 gsl::zstring<NI_MAXHOST> name) noexcept;

_INTERFACE_
std::errc nameof(const sockaddr_in6& ep, //
                 gsl::zstring<NI_MAXHOST> name) noexcept;

_INTERFACE_
std::errc nameof(const sockaddr_in6& ep, //
                 gsl::zstring<NI_MAXHOST> name,
                 gsl::zstring<NI_MAXSERV> serv) noexcept;

#endif // COROUTINE_NET_IO_H

ext: h

max_stars_repo_name: eddeighton/coroutine

max_issues_repo_name: eddeighton/coroutine

max_forks_repo_name: eddeighton/coroutine


text:
#include <petsc.h>
#include "compressibleFlow.h"
#include "mesh.h"

#define DIM 2                   /* Geometric dimension */

typedef struct {
    PetscReal gamma;
    PetscReal rhoL;
    PetscReal rhoR;
    PetscReal uL;
    PetscReal uR;
    PetscReal pL;
    PetscReal pR;
    PetscReal maxTime;
    PetscReal length;
} InitialConditions;

typedef struct {
    PetscReal rho;
    PetscReal rhoU[2];//Hardcode for 2 for now
    PetscReal rhoE;
} EulerNode;

typedef struct {
    PetscReal pstar,ustar,rhostarL,astarL,SL,SHL,STL,rhostarR,astarR,SR,SHR,STR,gamm1,gamp1;
} StarState;

typedef struct {
    InitialConditions initialConditions;
    StarState starState;
} ProblemSetup;

typedef struct{
    PetscReal f;
    PetscReal fprm;
} PressureFunction;

static PressureFunction f_and_fprm_rarefaction(PetscReal pstar, PetscReal pLR, PetscReal aLR, PetscReal gam, PetscReal gamm1,PetscReal gamp1) {
    // compute value of pressure function for rarefaction
    PressureFunction function;
    function.f = ((2. * aLR) / gamm1) * (pow(pstar / pLR, 0.5 * gamm1 / gam) - 1.);
    function.fprm = (aLR / pLR / gam) * pow(pstar / pLR, -0.5 * gamp1 / gam);
    return function;
}

static PressureFunction f_and_fprm_shock(PetscReal pstar, PetscReal pLR, PetscReal rhoLR, PetscReal gam, PetscReal gamm1, PetscReal gamp1){
    // compute value of pressure function for shock
    PetscReal A = 2./gamp1/rhoLR;
    PetscReal B = gamm1*pLR/gamp1;
    PetscReal sqrtterm =  PetscSqrtReal(A/(pstar+B));
    PressureFunction function;
    function.f = (pstar-pLR)*sqrtterm;
    function.fprm = sqrtterm*(1.-0.5*(pstar-pLR)/(B+pstar));
    return function;
}

#define EPS 1.e-6
#define MAXIT 100
PetscErrorCode DetermineStarState(const InitialConditions* setup, StarState* starState){
    // compute the speed of sound
    PetscReal aL = PetscSqrtReal(setup->gamma*setup->pL/setup->rhoL);
    PetscReal aR = PetscSqrtReal(setup->gamma*setup->pR/setup->rhoR);

    //first guess pstar based on two-rarefacation approximation
    starState->pstar = aL+aR - 0.5*(setup->gamma-1.)*(setup->uR-setup->uL);
    starState->pstar  = starState->pstar  / (aL/pow(setup->pL,0.5*(setup->gamma-1.)/setup->gamma) + aR/pow(setup->pR,0.5*(setup->gamma-1.)/setup->gamma) );
    starState->pstar = pow(starState->pstar,2.*setup->gamma/(setup->gamma-1.));
    starState->gamm1 = setup->gamma-1.;
    starState->gamp1 = setup->gamma+1.;

    PressureFunction fL;
    if (starState->pstar <= setup->pL){
        fL = f_and_fprm_rarefaction(starState->pstar, setup->pL,aL,setup->gamma,starState->gamm1,starState->gamp1);
    }else{
        fL = f_and_fprm_shock(starState->pstar,setup->pL,setup->rhoL,setup->gamma,starState->gamm1,starState->gamp1);
    }

    PressureFunction fR;
    if (starState->pstar <= setup->pR) {
        fR = f_and_fprm_rarefaction(starState->pstar, setup->pR, aR, setup->gamma, starState->gamm1, starState->gamp1);
    }else {
        fR = f_and_fprm_shock(starState->pstar, setup->pR, setup->rhoR, setup->gamma, starState->gamm1, starState->gamp1);
    }
    PetscReal delu = setup->uR-setup->uL;

    // iterate using Newton-Rapson
    if ((fL.f+fR.f+delu)> EPS) {
        // iterate using Newton-Rapson
        for(PetscInt it =0; it < MAXIT+4; it++){
            PetscReal pold = starState->pstar;
            starState->pstar = pold - (fL.f+fR.f+delu)/(fL.fprm+fR.fprm);

            if(starState->pstar < 0){
                starState->pstar = EPS;
            }

            if(2.0*PetscAbsReal(starState->pstar - pold)/(starState->pstar + pold) < EPS){
                break;
            }else{
                if(starState->pstar < setup->pL){
                    fL = f_and_fprm_rarefaction(starState->pstar, setup->pL,aL,setup->gamma,starState->gamm1,starState->gamp1);
                }else{
                    fL = f_and_fprm_shock(starState->pstar,setup->pL,setup->rhoL,setup->gamma,starState->gamm1,starState->gamp1);
                }
                if (starState->pstar<=setup->pR) {
                    fR = f_and_fprm_rarefaction(starState->pstar, setup->pR, aR, setup->gamma, starState->gamm1, starState->gamp1);
                }else {
                    fR = f_and_fprm_shock(starState->pstar, setup->pR, setup->rhoR, setup->gamma, starState->gamm1, starState->gamp1);
                }
            }

            if (it>MAXIT){
                SETERRQ(PETSC_COMM_WORLD,1,"error in Riemann.find_pstar - did not converage for pstar" );
            }
        }
    }

    // determine rest of star state
    starState->ustar = 0.5*(setup->uL+setup->uR+fR.f-fL.f);

    // left star state
    PetscReal pratio = starState->pstar/setup->pL;
    if (starState->pstar<=setup->pL) {  // rarefaction
        starState->rhostarL = setup->rhoL * PetscPowReal(pratio, 1. / setup->gamma);
        starState->astarL = aL * PetscPowReal(pratio, 0.5 * starState->gamm1 / setup->gamma);
        starState->SHL = setup->uL - aL;
        starState->STL = starState->ustar - starState->astarL;
    }else {  // #shock
        starState->rhostarL = setup->rhoL * (pratio + starState->gamm1 / starState->gamp1) / (starState->gamm1 * pratio / starState->gamp1 + 1.);
        starState->SL = setup->uL - aL * PetscSqrtReal(0.5 * starState->gamp1 / setup->gamma * pratio + 0.5 * starState->gamm1 / setup->gamma);
    }

    // right star state
    pratio = starState->pstar/setup->pR;
    if (starState->pstar<=setup->pR) {  // # rarefaction
        starState->rhostarR = setup->rhoR * PetscPowReal(pratio, 1. / setup->gamma);
        starState->astarR = aR * PetscPowReal(pratio, 0.5 * starState->gamm1 / setup->gamma);
        starState-> SHR = setup->uR + aR;
        starState->STR = starState->ustar + starState->astarR;
    }else {  // shock
        starState->rhostarR = setup->rhoR * (pratio + starState->gamm1 / starState->gamp1) / (starState->gamm1 * pratio / starState->gamp1 + 1.);
        starState->SR = setup->uR + aR * PetscSqrtReal(0.5 * starState->gamp1 / setup->gamma * pratio + 0.5 * starState->gamm1 / setup->gamma);
    }
    return 0;
}

void SetExactSolutionAtPoint(PetscInt dim, PetscReal xDt, const InitialConditions* setup, const StarState* starState, EulerNode* uu){
    PetscReal p;
    // compute the speed of sound
    PetscReal aL = PetscSqrtReal(setup->gamma*setup->pL/setup->rhoL);
    PetscReal aR = PetscSqrtReal(setup->gamma*setup->pR/setup->rhoR);

    for(PetscInt i =0; i < dim; i++){
        uu->rhoU[i] = 0.0;
    }

    if (xDt <= starState->ustar) {  //# left of contact surface
        if (starState->pstar <= setup->pL) {  // # rarefaction
            if (xDt <= starState->SHL) {
                uu->rho = setup->rhoL;
                p = setup->pL;
                uu->rhoU[0] = setup->uL*uu->rho;
            }else if (xDt <=starState->STL) {  //#SHL < x / t < STL
                PetscReal tmp = 2. / starState->gamp1 + (starState->gamm1 / starState->gamp1 / aL) * (setup->uL - xDt);
                uu->rho = setup->rhoL * pow(tmp, 2. / starState->gamm1);
                uu->rhoU[0] = uu->rho * (2. / starState->gamp1) * (aL + 0.5 * starState->gamm1 * setup->uL + xDt);
                p = setup->pL * pow(tmp, 2. * setup->gamma / starState->gamm1);
            }else {  //# STL < x/t < u*
                uu->rho = starState->rhostarL;
                p = starState->pstar;
                uu->rhoU[0] = uu->rho * starState->ustar;
            }
        }else{  //# shock
            if (xDt<= starState->SL) {  // # xDt < SL
                uu->rho = setup->rhoL;
                p = setup->pL;
                uu->rhoU[0] = uu->rho * setup->uL;
            }else {  //# SL < xDt < ustar
                uu->rho = starState->rhostarL;
                p = starState->pstar;
                uu->rhoU[0] = uu->rho * starState->ustar;
            }
        }
    }else{//# right of contact surface
        if (starState->pstar<=setup->pR) {  //# rarefaction
            if (xDt>= starState->SHR) {
                uu->rho = setup->rhoR;
                p = setup->pR;
                uu->rhoU[0] = uu->rho * setup->uR;
            }else if (xDt >= starState->STR) {  // # SHR < x/t < SHR
                PetscReal tmp = 2./starState->gamp1 - (starState->gamm1/starState->gamp1/aR)*(setup->uR-xDt);
                uu->rho = setup->rhoR*PetscPowReal(tmp,2./starState->gamm1);
                uu->rhoU[0] = uu->rho * (2./starState->gamp1)*(-aR + 0.5*starState->gamm1*setup->uR+xDt);
                p = setup->pR*PetscPowReal(tmp,2.*setup->gamma/starState->gamm1);
            }else{ //# u* < x/t < STR
                uu->rho = starState->rhostarR;
                p = starState->pstar;
                uu->rhoU[0] = uu->rho * starState->ustar;
            }
        }else {//# shock
            if (xDt>= starState->SR) {  // # xDt > SR
                uu->rho = setup->rhoR;
                p = setup->pR;
                uu->rhoU[0] = uu->rho * setup->uR;
            }else {//#ustar < xDt < SR
                uu->rho = starState->rhostarR;
                p = starState->pstar;
                uu->rhoU[0] = uu->rho * starState->ustar;
            }
        }
    }
    PetscReal e =  p/starState->gamm1/uu->rho;
    PetscReal E = e + 0.5*(uu->rhoU[0]/uu->rho)*(uu->rhoU[0]/uu->rho);
    uu->rhoE = uu->rho*E;
}

static PetscErrorCode SetExactSolutionRho(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx){
    ProblemSetup* prob = (ProblemSetup*)ctx;

    PetscReal xDt = (x[0]-prob->initialConditions.length/2)/time;
    EulerNode uu;
    SetExactSolutionAtPoint(dim, xDt, &prob->initialConditions, &prob->starState, &uu);

    u[0] = uu.rho;
    return 0;
}

static PetscErrorCode SetExactSolutionRhoU(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx){
    ProblemSetup* prob = (ProblemSetup*)ctx;

    PetscReal xDt = (x[0]-prob->initialConditions.length/2)/time;
    EulerNode uu;
    SetExactSolutionAtPoint(dim, xDt, &prob->initialConditions, &prob->starState, &uu);
    u[0] = uu.rhoU[0];
    u[1] = uu.rhoU[1];
    return 0;
}

static PetscErrorCode SetExactSolutionRhoE(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx){
    ProblemSetup* prob = (ProblemSetup*)ctx;

    PetscReal xDt = (x[0]-prob->initialConditions.length/2)/time;
    EulerNode uu;
    SetExactSolutionAtPoint(dim, xDt, &prob->initialConditions, &prob->starState, &uu);
    u[0] = uu.rhoE;
    return 0;
}

static PetscErrorCode PrintVector(DM dm, Vec v, PetscInt step, const char * fileName){
    Vec                cellgeom;
    PetscErrorCode ierr = DMPlexGetGeometryFVM(dm, NULL, &cellgeom, NULL);CHKERRQ(ierr);
    PetscInt cStart, cEnd;
        ierr = DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
    DM dmCell;
    ierr = VecGetDM(cellgeom, &dmCell);CHKERRQ(ierr);
    const PetscScalar *cgeom;
    ierr = VecGetArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
    const PetscScalar      *x;
    ierr = VecGetArrayRead(v, &x);CHKERRQ(ierr);
    // print the header for each file
    char filename[100];
    ierr = PetscSNPrintf(filename,sizeof(filename),"%s.%d.txt",fileName, step);CHKERRQ(ierr);

    PetscInt rank = 0;
    PetscInt size;
    ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRMPI(ierr);
    ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);CHKERRMPI(ierr);

    for(PetscInt r =0; r < size; r++ ) {
        if(r == rank) {
            FILE *fptr;
            if(r == 0){
                fptr = fopen(filename, "w");
                fprintf(fptr, "x rho u e rhou rhoe\n");
            }else{
                fptr = fopen(filename, "a");
            }
            for (PetscInt c = cStart; c < cEnd; ++c) {
                PetscFVCellGeom *cg;
                const EulerNode *xc;

                ierr = DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
                CHKERRQ(ierr);
                ierr = DMPlexPointGlobalFieldRead(dm, c, 0, x, &xc);
                CHKERRQ(ierr);
                if(xc) {// must be real cell and not ghost
                    PetscReal u0 = xc->rhoU[0] / xc->rho;
                    fprintf(fptr, "%f %f %f %f %f %f\n", cg->centroid[0], xc->rho, u0, (xc->rhoE / xc->rho) - 0.5 * u0 * u0, xc->rhoU[0], xc->rhoE);
                }
            }

            fclose(fptr);
        }
        MPI_Barrier(PetscObjectComm((PetscObject)dm));
    }
    ierr = VecRestoreArrayRead(cellgeom, &cgeom);CHKERRQ(ierr);
    ierr = VecRestoreArrayRead(v, &x);CHKERRQ(ierr);
    return 0;
}

static PetscErrorCode MonitorError(TS ts, PetscInt step, PetscReal time, Vec u, void *ctx) {
    PetscFunctionBeginUser;
    PetscErrorCode     ierr;

    // Get the DM
    DM dm;
    ierr = TSGetDM(ts, &dm);CHKERRQ(ierr);

    // Create a copy of the u vector
    Vec e;
    ierr = DMCreateGlobalVector(dm, &e);CHKERRQ(ierr);
    ierr = PetscObjectSetName((PetscObject)e, "exact");CHKERRQ(ierr);

    // Set the values
    PetscErrorCode     (*func[3]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {SetExactSolutionRho, SetExactSolutionRhoU, SetExactSolutionRhoE};
    void* ctxs[3] ={ctx, ctx, ctx};
    ierr    = DMProjectFunction(dm,time,func,ctxs,INSERT_ALL_VALUES,e);CHKERRQ(ierr);

    // just print to a file for now
    ierr = PrintVector(dm, e, step, "exact");CHKERRQ(ierr);
    ierr = PrintVector(dm, u, step, "solution");CHKERRQ(ierr);

    PetscPrintf(PETSC_COMM_WORLD, "TS %d: %f\n", step, time);

    DMRestoreGlobalVector(dm, &e);
    PetscFunctionReturn(0);
}

static PetscErrorCode PhysicsBoundary_Euler_Mirror(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)
{
    const EulerNode *xI = (const EulerNode*)a_xI;
    EulerNode       *xG = (EulerNode*)a_xG;
    ProblemSetup* prob = (ProblemSetup*)ctx;
    PetscFunctionBeginUser;
    xG->rho = xI->rho;
    xG->rhoE = xI->rhoE;
    xG->rhoU[0] = xI->rhoU[0];
    xG->rhoU[1] = xI->rhoU[1];

    PetscFunctionReturn(0);
}

/* PetscReal* => EulerNode* conversion */
static PetscErrorCode PhysicsBoundary_Euler_Left(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)
{
    const EulerNode *xI = (const EulerNode*)a_xI;
    EulerNode       *xG = (EulerNode*)a_xG;
    ProblemSetup* prob = (ProblemSetup*)ctx;
    PetscFunctionBeginUser;
    xG->rho = prob->initialConditions.rhoL;
    PetscReal eT = prob->initialConditions.rhoL*((prob->initialConditions.pL /(prob->initialConditions.gamma -1) / prob->initialConditions.rhoL) + 0.5 * prob->initialConditions.uL * prob->initialConditions.uL);
    xG->rhoE = eT;
    xG->rhoU[0] = prob->initialConditions.rhoL * prob->initialConditions.uL;
    xG->rhoU[1] = 0.0;

    PetscFunctionReturn(0);
}

/* PetscReal* => EulerNode* conversion */
static PetscErrorCode PhysicsBoundary_Euler_Right(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)
{
    const EulerNode *xI = (const EulerNode*)a_xI;
    EulerNode       *xG = (EulerNode*)a_xG;
    ProblemSetup*        prob = (ProblemSetup*)ctx;
    PetscFunctionBeginUser;
    xG->rho = prob->initialConditions.rhoR;
    PetscReal eT = prob->initialConditions.rhoR*((prob->initialConditions.pR /(prob->initialConditions.gamma -1)/ prob->initialConditions.rhoR) + 0.5 * prob->initialConditions.uR * prob->initialConditions.uR);
    xG->rhoE = eT;
    xG->rhoU[0] = prob->initialConditions.rhoR * prob->initialConditions.uR;
    xG->rhoU[1] = 0.0;

    PetscFunctionReturn(0);
}

int main(int argc, char **argv)
{
    PetscErrorCode ierr;
    // create the mesh
    // setup the ts
    DM dm;                 /* problem definition */
    TS ts;                 /* timestepper */

    // initialize petsc and mpi
    PetscInitialize(&argc, &argv, NULL, "HELP");

    // Create a ts
    ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr);
    ierr = TSSetProblemType(ts, TS_NONLINEAR);CHKERRQ(ierr);
    ierr = TSSetType(ts, TSEULER);CHKERRQ(ierr);
    ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);

    // Create a mesh
    // hard code the problem setup
    PetscReal start[] = {0.0, 0.0};
    PetscReal end[] = {1.0, 1};
    PetscInt nx[] = {100, 1};
    DMBoundaryType bcType[] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
    ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, DIM, PETSC_FALSE, nx, start, end, bcType, PETSC_TRUE, &dm);CHKERRQ(ierr);

// Setup the problem
    ProblemSetup problem;

    // case 1 - Sod problem
    problem.initialConditions.rhoL=1.0;
    problem.initialConditions.uL=0.0;
    problem.initialConditions.pL=1.0;
    problem.initialConditions.rhoR=0.125;
    problem.initialConditions.uR=0.0;
    problem.initialConditions.pR=0.1;
    problem.initialConditions.maxTime = 0.25;
    problem.initialConditions.length = 1;
    problem.initialConditions.gamma = 1.4;

    // case 2
//    problem.initialConditions.rhoL=1.0;
//    problem.initialConditions.uL=-2.0;
//    problem.initialConditions.pL=0.4;
//    problem.initialConditions.rhoR=1.0;
//    problem.initialConditions.uR=2.0;
//    problem.initialConditions.pR=0.4;
//    problem.initialConditions.maxTime = 0.15;
//    problem.initialConditions.length = 1;
//    problem.initialConditions.gamma = 1.4;

    // case 5
//    problem.initialConditions.rhoL=5.99924;
//    problem.initialConditions.uL=19.5975;
//    problem.initialConditions.pL=460.894;
//    problem.initialConditions.rhoR=5.99242;
//    problem.initialConditions.uR=-6.19633;
//    problem.initialConditions.pR=46.0950;
//    problem.initialConditions.maxTime =  0.035;
//    problem.initialConditions.length = 1;
//    problem.initialConditions.gamma = 1.4;

    // Compute the star state
    ierr = DetermineStarState(&problem.initialConditions, &problem.starState);CHKERRQ(ierr);

    // Setup the flow data
    FlowData flowData;     /* store some of the flow data*/
    ierr = FlowCreate(&flowData);CHKERRQ(ierr);

    //Setup
    CompressibleFlow_SetupDiscretization(flowData, dm);

    // Add in the flow parameters
    PetscScalar params[TOTAL_COMPRESSIBLE_FLOW_PARAMETERS];
    params[CFL] = 0.5;
    params[GAMMA] = problem.initialConditions.gamma;

    // set up the finite volume fluxes
    CompressibleFlow_StartProblemSetup(flowData, TOTAL_COMPRESSIBLE_FLOW_PARAMETERS, params);
    DMView(flowData->dm, PETSC_VIEWER_STDERR_SELF);
    // Add in any boundary conditions
    PetscDS prob;
    ierr = DMGetDS(flowData->dm, &prob);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    const PetscInt idsLeft[]= {4};
    ierr = PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall left", "Face Sets", 0, 0, NULL, (void (*)(void)) PhysicsBoundary_Euler_Left, NULL, 1, idsLeft, &problem);CHKERRQ(ierr);
    const PetscInt idsRight[]= {2};
    ierr = PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall right", "Face Sets", 0, 0, NULL, (void (*)(void)) PhysicsBoundary_Euler_Right, NULL, 1, idsRight, &problem);CHKERRQ(ierr);
    const PetscInt mirror[]= {1, 3};
    ierr = PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "mirrorWall", "Face Sets", 0, 0, NULL, (void (*)(void)) PhysicsBoundary_Euler_Mirror, NULL, 2, mirror, &problem);CHKERRQ(ierr);

    // Complete the problem setup
    ierr = CompressibleFlow_CompleteProblemSetup(flowData, ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Name the flow field
    ierr = PetscObjectSetName(((PetscObject)flowData->flowField), "Numerical Solution");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Setup the TS
    ierr = TSSetFromOptions(ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSMonitorSet(ts, MonitorError, &problem, NULL);CHKERRQ(ierr);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSSetMaxTime(ts,problem.initialConditions.maxTime);CHKERRQ(ierr);

    // set the initial conditions
    PetscErrorCode     (*func[3]) (PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {SetExactSolutionRho, SetExactSolutionRhoU, SetExactSolutionRhoE};
    void* ctxs[3] ={&problem, &problem, &problem};
    ierr    = DMProjectFunction(flowData->dm,0.0,func,ctxs,INSERT_ALL_VALUES,flowData->flowField);CHKERRQ(ierr);

    ierr = TSSolve(ts,flowData->flowField);CHKERRQ(ierr);

    return PetscFinalize();

}

ext: c

max_stars_repo_name: mmcgurn/MattFlowCases

max_issues_repo_name: mmcgurn/MattFlowCases

max_forks_repo_name: mmcgurn/MattFlowCases


text:
/**
 *
 * @file testing_zpotri.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_zmain.h"

enum blas_order_type {
            blas_rowmajor = 101,
            blas_colmajor = 102 };

enum blas_cmach_type {
            blas_base      = 151,
            blas_t         = 152,
            blas_rnd       = 153,
            blas_ieee      = 154,
            blas_emin      = 155,
            blas_emax      = 156,
            blas_eps       = 157,
            blas_prec      = 158,
            blas_underflow = 159,
            blas_overflow  = 160,
            blas_sfmin     = 161};

enum blas_norm_type {
            blas_one_norm       = 171,
            blas_real_one_norm  = 172,
            blas_two_norm       = 173,
            blas_frobenius_norm = 174,
            blas_inf_norm       = 175,
            blas_real_inf_norm  = 176,
            blas_max_norm       = 177,
            blas_real_max_norm  = 178 };

static void
BLAS_error(char *rname, int err, int val, int x) {
  fprintf( stderr, "%s %d %d %d\n", rname, err, val, x );
  abort();
}

static
void
BLAS_zge_norm(enum blas_order_type order, enum blas_norm_type norm,
  int m, int n, const PLASMA_Complex64_t *a, int lda, double *res) {
  int i, j; float anorm, v;
  char rname[] = "BLAS_zge_norm";

  if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 );

  if (norm == blas_frobenius_norm) {
    anorm = 0.0f;
    for (j = n; j; --j) {
      for (i = m; i; --i) {
        v = a[0];
        anorm += v * v;
        a++;
      }
      a += lda - m;
    }
    anorm = sqrt( anorm );
  } else if (norm == blas_inf_norm) {
    anorm = 0.0f;
    for (i = 0; i < m; ++i) {
      v = 0.0f;
      for (j = 0; j < n; ++j) {
        v += cabs( a[i + j * lda] );
      }
      if (v > anorm)
        anorm = v;
    }
  } else if (norm == blas_one_norm) {
    anorm = 0.0f;
    for (i = 0; i < m; ++i) {
      v = 0.0f;
      for (j = 0; j < n; ++j) {
        v += cabs( a[i + j * lda] );
      }
      if (v > anorm)
        anorm = v;
    }
  } else {
    BLAS_error( rname, -2, norm, 0 );
    return;
  }

  if (res) *res = anorm;
}

static
double
BLAS_dpow_di(double x, int n) {
  double rv = 1.0;

  if (n < 0) {
    n = -n;
    x = 1.0 / x;
  }

  for (; n; n >>= 1, x *= x) {
    if (n & 1)
      rv *= x;
  }

  return rv;
}

static
double
BLAS_dfpinfo(enum blas_cmach_type cmach) {
  double eps = 1.0, r = 1.0, o = 1.0, b = 2.0;
  int t = 53, l = 1024, m = -1021;
  char rname[] = "BLAS_dfpinfo";

  if ((sizeof eps) == sizeof(float)) {
    t = 24;
    l = 128;
    m = -125;
  } else {
    t = 53;
    l = 1024;
    m = -1021;
  }

  /* for (i = 0; i < t; ++i) eps *= half; */
  eps = BLAS_dpow_di( b, -t );
  /* for (i = 0; i >= m; --i) r *= half; */
  r = BLAS_dpow_di( b, m-1 );

  o -= eps;
  /* for (i = 0; i < l; ++i) o *= b; */
  o = (o * BLAS_dpow_di( b, l-1 )) * b;

  switch (cmach) {
    case blas_eps: return eps;
    case blas_sfmin: return r;
    default:
      BLAS_error( rname, -1, cmach, 0 );
      break;
  }
  return 0.0;
}

static int check_factorization(int, PLASMA_Complex64_t*, PLASMA_Complex64_t*, int, int, double);
static int check_inverse(int, PLASMA_Complex64_t *, PLASMA_Complex64_t *, int, int, double);

int testing_zpotri(int argc, char **argv)
{

    /* Check for number of arguments*/
    if (argc != 2){
        USAGE("POTRI", "N LDA",
              "   - N    : the size of the matrix\n"
              "   - LDA  : leading dimension of the matrix A\n");
        return -1;
    }

    int N     = atoi(argv[0]);
    int LDA   = atoi(argv[1]);
    double eps;
    int uplo;
    int info_inverse, info_factorization;

    PLASMA_Complex64_t *A1   = (PLASMA_Complex64_t *)malloc(LDA*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *A2   = (PLASMA_Complex64_t *)malloc(LDA*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *WORK = (PLASMA_Complex64_t *)malloc(2*LDA*sizeof(PLASMA_Complex64_t));

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)){
        printf("Out of Memory \n ");
        return -2;
    }

    eps = BLAS_dfpinfo( blas_eps );

    uplo = PlasmaUpper;

    /*-------------------------------------------------------------
    *  TESTING ZPOTRI
    */

    /* Initialize A1 and A2 for Symmetric Positif Matrix */
    PLASMA_zplghe( (double)N, N, A1, LDA, 51 );
    PLASMA_zlacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

    printf("\n");
    printf("------ TESTS FOR PLASMA ZPOTRI ROUTINE -------  \n");
    printf("            Size of the Matrix %d by %d\n", N, N);
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n", eps);
    printf(" Computational tests pass if scaled residuals are less than 60.\n");

    /* PLASMA ZPOTRF */
    PLASMA_zpotrf(uplo, N, A2, LDA);

    /* Check the factorization */
    info_factorization = check_factorization( N, A1, A2, LDA, uplo, eps);

    /* PLASMA ZPOTRI */
    PLASMA_zpotri(uplo, N, A2, LDA);

    /* Check the inverse */
    info_inverse = check_inverse(N, A1, A2, LDA, uplo, eps);

    if ( (info_inverse == 0) && (info_factorization == 0) ) {
        printf("***************************************************\n");
        printf(" ---- TESTING ZPOTRI ..................... PASSED !\n");
        printf("***************************************************\n");
    }
    else {
        printf("***************************************************\n");
        printf(" - TESTING ZPOTRI ... FAILED !\n");
        printf("***************************************************\n");
    }

    free(A1); free(A2); free(WORK);

    return 0;
}

/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */
static int check_factorization(int N, PLASMA_Complex64_t *A1, PLASMA_Complex64_t *A2, int LDA, int uplo, double eps)
{
    double Anorm, Rnorm;
    PLASMA_Complex64_t alpha;
    int info_factorization;
    int i,j;

    PLASMA_Complex64_t *Residual = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *L1       = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *L2       = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    double *work              = (double *)malloc(N*sizeof(double));

    memset((void*)L1, 0, N*N*sizeof(PLASMA_Complex64_t));
    memset((void*)L2, 0, N*N*sizeof(PLASMA_Complex64_t));

    alpha= 1.0;

    LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,' ', N, N, A1, LDA, Residual, N);

    /* Dealing with L'L or U'U  */
    if (uplo == PlasmaUpper){
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L1, N);
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L2, N);
        cblas_ztrmm(CblasColMajor, CblasLeft, CblasUpper, CblasConjTrans, CblasNonUnit, N, N, CBLAS_SADDR(alpha), L1, N, L2, N);
    }
    else{
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L1, N);
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L2, N);
        cblas_ztrmm(CblasColMajor, CblasRight, CblasLower, CblasConjTrans, CblasNonUnit, N, N, CBLAS_SADDR(alpha), L1, N, L2, N);
    }

    /* Compute the Residual || A -L'L|| */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
           Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];

    BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, N, Residual, N, &Rnorm );
    BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, N, A1, LDA, &Anorm );

    printf("============\n");
    printf("Checking the Cholesky Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));

    if ( isnan(Rnorm/(Anorm*N*eps)) || isinf(Rnorm/(Anorm*N*eps)) || (Rnorm/(Anorm*N*eps) > 60.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}

/*------------------------------------------------------------------------
 *  Check the accuracy of the computed inverse
 */

static int check_inverse(int N, PLASMA_Complex64_t *A1, PLASMA_Complex64_t *A2, int LDA, int uplo, double eps )
{
    int info_inverse;
    int i, j;
    double Rnorm, Anorm, Ainvnorm, result;
    PLASMA_Complex64_t alpha, beta, zone;
    PLASMA_Complex64_t *work = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));

    alpha = -1.0;
    beta  = 0.0;
    zone = 1.0;

    /* Rebuild the other part of the inverse matrix */
    if(uplo == PlasmaUpper){
       for(j=0; j<N; j++)
          for(i=0; i<j; i++)
             *(A2+j+i*LDA) = *(A2+i+j*LDA);
       cblas_zhemm(CblasColMajor, CblasLeft, CblasUpper, N, N, CBLAS_SADDR(alpha), A2, LDA, A1, LDA, CBLAS_SADDR(beta), work, N);

    }
    else {
       for(j=0; j<N; j++)
          for(i=j; i<N; i++)
             *(A2+j+i*LDA) = *(A2+i+j*LDA);
       cblas_zhemm(CblasColMajor, CblasLeft, CblasLower, N, N, CBLAS_SADDR(alpha), A2, LDA, A1, LDA, CBLAS_SADDR(beta), work, N);
    }

    /* Add the identity matrix to work */
    for(i=0; i<N; i++)
        *(work+i+i*N) = *(work+i+i*N) + zone;

    BLAS_zge_norm( blas_colmajor, blas_one_norm, N, N, work, N, &Rnorm );
    BLAS_zge_norm( blas_colmajor, blas_one_norm, N, N, A1, LDA, &Anorm );
    BLAS_zge_norm( blas_colmajor, blas_one_norm, N, N, A2, LDA, &Ainvnorm );

    if (getenv("PLASMA_TESTING_VERBOSE"))
      printf( "||A||_1=%f\n||Ainv||_1=%f\n||Id - A*Ainv||_1=%e\n", Anorm, Ainvnorm, Rnorm );

    result = Rnorm / ( (Anorm*Ainvnorm)*N*eps ) ;
    printf("============\n");
    printf("Checking the Residual of the inverse \n");
    printf("-- ||Id - A*Ainv||_1/((||A||_1||Ainv||_1).N.eps) = %e \n", result);

    if (  isnan(Ainvnorm) || isinf(Ainvnorm) || isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- The inverse is suspicious ! \n");
        info_inverse = 1;
     }
    else{
        printf("-- The inverse is CORRECT ! \n");
        info_inverse = 0;
    }

    free(work);

    return info_inverse;
}