#define RLC_Z0  50.0f
#define RLC_PI  3.14159265358979323846f
#define RLC_X_MIN   1.0f

/* globali già definite altrove:
   float vna_s11_re[VNA_MAX_POINTS];
   float vna_s11_im[VNA_MAX_POINTS];
   U32   vna_frequency[VNA_MAX_POINTS];
   U16   vna_points;
   char  rlc;       // 0=R, 1=L, 2=C
   float rlc_val;   // ohm, nH, pF
*/

#include <math.h>

// =====================================================
// cerca:
// rlc  0=r 1=l 2=c
// best  miglior punto
// quale guale a best
// val valore (Ohm, pF, nHe)
// =====================================================

int do_rlc(void)
{
    int i;
    int j;
    int start=-1;
    int stop=-1;
    int best=-1;
    int type=-1;
    float best_err=1.0e30f;
    float gr,gi;
    float den;
    float r,x;
    float val;
    double f;

    rlc_val=0.0f;
    if(vna_points < 5) return 0;

// =====================================================
// 1) identifica il tipo nella prima zona utile
//    a bassa frequenza
// =====================================================
//  RESISTENZA ?
// =====================================================
//  valor resistivo a basse frequenze
    rlc_val=0;
    float rlc_z=0;
    for(int j=1;j<4;j++)
    {
        rlc_val += vna_z_re[j];
        rlc_z   += vna_z_im[j];
    }
    rlc_val /= 3.;
    rlc_z   /= 3.;
    if(rlc_val > fabsf(rlc_z))
    if (rlc_val > 0.9)
    {
       int n=0;
       rlc_val=0;
       for(int j=1;j<11 && j<vna_points;j++)
       {
          rlc_val += vna_z_re[j];
          n++;
       }
       if(n==0) return 0;
       rlc_val /= (float)n;
       rlc=0;
       type=0;
       start=1;
       quale=1;
       best=1;
       if (rlc_val < 0) return 0;
       if (fabsf(rlc_val) > 20000.0f) return 0;
       freq_best=vna_frequency[best]/1000000.;
       return 1;
    }
// =====================================================
// Casi INDUTTANZA e CAPACITA'
// =====================================================
    for(i=0;i<vna_points;i++)
    {
        gr=vna_s11_re[i];
        gi=vna_s11_im[i];
        den=(1.0f-gr)*(1.0f-gr)+gi*gi;
//      troppo vicino ad OPEN
        if(den < 1.0e-8f) continue;
//      Z = 50 * (1+Gamma)/(1-Gamma)
        r=RLC_Z0 * (1.0f-gr*gr-gi*gi) / den;
        x=RLC_Z0 * (2.0f*gi) / den;
        if(!isfinite(r) || !isfinite(x)) continue;
//      Evita valori estremi dove la misura
//      è poco significativa.
        if(fabsf(x) > 5000.0f) continue;
// =====================================================
//      INDUTTANZA ?
// =====================================================
        if(x > 0.5f)
        {
            type=1;
            start=i;
            break;
        }
// =====================================================
//      CAPACITA' ?
// =====================================================
        if(x>=0.0f) continue;
        if(fabsf(gi) < 0.01f) continue;
        if(x < -0.5f)
        {
            type=2;
            start=i;
            break;
        }
    }
    if(start < 0)
    return 0;
// =====================================================
// 2) trova il limite della prima zona utile
//    prima della prima risonanza
// =====================================================
    stop=vna_points-1;
    if(type==1 || type==2)
    {
        for(i=start+1;i<vna_points;i++)
        {
            gr=vna_s11_re[i];
            gi=vna_s11_im[i];
            den=(1.0f-gr)*(1.0f-gr)+gi*gi;
            if(den < 1.0e-8f) continue;
            r=RLC_Z0 * (1.0f-gr*gr-gi*gi) / den;
            x=RLC_Z0 * (2.0f*gi) / den;
            if(!isfinite(r) || !isfinite(x)) continue;
//          Prima risonanza di un induttore:
//          X passa da positivo a zero/negativo.
            if(type==1 && x<=0.0f)
            {
                stop=i-1;
                break;
            }
//          Prima risonanza di un condensatore:
//          X passa da negativo a zero/positivo.
            if(type==2 && x>=0.0f)
            {
                stop=i-1;
                break;
            }
        }
    }
    if(stop < start) return 0;
//  =====================================================
//  3) cerca il punto migliore
//  =====================================================
    for(int j=start;j<=stop;j++)
    {
//      if(type==0) i=start+4; else
        i=j;
        float err;
        gr=vna_s11_re[i];
        gi=vna_s11_im[i];
        den=(1.0f-gr)*(1.0f-gr)+gi*gi;
        if(den < 1.0e-8f) continue;
        r=RLC_Z0 * (1.0f-gr*gr-gi*gi) / den;
        x=RLC_Z0 * (2.0f*gi) / den;
        if(!isfinite(r) || !isfinite(x)) continue;
        f=(double)vna_frequency[i];
        if(f<=0.0) continue;
// -------------------------------------------------
//     INDUTTANZA
// -------------------------------------------------
        if(type==1)
        {
            if(x<=0.0f)  continue;
//          L in nH
            val=(float) (((double)x / (2.0*RLC_PI*f)) *1.0e9 );
//          Preferisce X vicino a 50 ohm.
            err=fabsf(fabsf(x)-50.0f);
        }
// -------------------------------------------------
//      CAPACITA'
// -------------------------------------------------
        else if(type==2)
        {
            if(x>=0.0f) continue;
//          C in pF
            val=(float) ((-1.0 /(2.0*RLC_PI*f*(double)x))*1.0e12);
//          Preferisce |X| vicino a 50 ohm.
            err=fabsf(fabsf(x)-50.0f);
        }
        if(!isfinite(val) || val<=0.0f) continue;
        if(err < best_err)
        {
            best_err=err;
            best=i;
            quale=i;
            rlc_val=val;
        }
    }
// =====================================================
    if(best < 0) return 0;
    if(best > 5)best--;
    if(best > 5)best--;
    if(best < start) best=start;
    quale=best;
    rlc=(char)type;
//  ricalcola il valore sul nuovo best
    gr=vna_s11_re[best];
    gi=vna_s11_im[best];
    den=(1.0f-gr)*(1.0f-gr)+gi*gi;
    if(den < 1.0e-8f) return 0;
    r=RLC_Z0 * (1.0f-gr*gr-gi*gi) / den;
    x=RLC_Z0 * (2.0f*gi) / den;
    if(!isfinite(r) || !isfinite(x)) return 0;
    f=(double)vna_frequency[best];
    if(f<=0.0) return 0;
    if(type==1)
    {
        if(x<=0.0f) return 0;
        rlc_val=(float) (((double)x /(2.0*RLC_PI*f))*1.0e9);
    }
    else if(type==2)
    {
        if(x>=0.0f) return 0;
        rlc_val=(float) ((-1.0 /(2.0*RLC_PI*f*(double)x))*1.0e12);
    }
    if(!isfinite(rlc_val) || rlc_val<=0.0f) return 0;

// =====================================================
// 4) media locale sui 5 punti attorno al migliore
//    vengono accettati solo valori entro +/-20%
//    rispetto al valore del punto migliore
// =====================================================
    {
        int n=0;
        int j0;
        int j1;
        double sum=0.0;
        float ref_val=rlc_val;
        j0=best-2;
        j1=best+2;
        if(j0<start)
        j0=start;
        if(j1>stop)
        j1=stop;
        for(j=j0;j<=j1;j++)
        {
            gr=vna_s11_re[j];
            gi=vna_s11_im[j];
            den=(1.0f-gr)*(1.0f-gr)+gi*gi;
            if(den < 1.0e-8f) continue;
            r=RLC_Z0 *
            (1.0f-gr*gr-gi*gi) / den;
             x=RLC_Z0 *
            (2.0f*gi) / den;
            if(!isfinite(r) || !isfinite(x)) continue;
            f=(double)vna_frequency[j];
            if(f<=0.0) continue;
// ---------------- L ----------------
            if(type==1)
            {
                if(x<=0.0f) continue;
                val=(float) (((double)x /(2.0*RLC_PI*f))*1.0e9);
            }
// ---------------- C ----------------
            else if(type==2)
            {
                if(x>=0.0f) continue;
                val=(float) ((-1.0 /(2.0*RLC_PI*f*(double)x))*1.0e12);
            }
            if(!isfinite(val) || val<=0.0f) continue;
//          Scarta punti che differiscono
//          più del 20% dal valore centrale.
            if(fabsf(val-ref_val) >
            fabsf(ref_val)*0.20f) continue;
            sum+=(double)val;
            n++;
        }
//      Se almeno due punti sono coerenti,
//      usa la loro media.
        if(n>=2) rlc_val=(float)(sum/(double)n);
    }
//  cap.parassite
    if (type==1) rlc_val -= 0;   //6.7 nHe compensate in SHORT
    if (type==2) rlc_val -= 0.1; //0.1 pF punta
    if (rlc_val < 0.5) return 0;
    freq_best=vna_frequency[best]/1000000.;
    quale=best;
    return 1;
}