Garnet

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
namespace Scientific
{
    class Ferrite
    {
        #region [ 参数 ]
        double pi = Math.PI;
        double U0;   //真空磁导率    u0
        double E0;   //真空介电常数  E0
        double r;    //旋磁比        r
        double K;
        protected double C;    //真空光速      C
        public double Ms;   //饱和磁场强度  4πMs
        public double Er;   //相对介电常数  Er
        public double losstan;
        public double deltaH;
        public double Ueff; //相对磁导率    Ueff
        public double ku_c;//由内导体计算出所需的k/u
        #endregion
        #region [ 构造函数 ]
        public Ferrite()
        {
            U0 = 4 * Math.PI * 1e-7;
            E0 = 8.854e-12;
            r = 2.8;
            C = 1 / Math.Sqrt(U0 * E0);
            Ueff = 1.0D;
            losstan = 0D;
            deltaH = 0D;
            K = 2 * pi * 2e6 / C;//2GHz时的真空K0;k=2*pi*f/c;
            ku_c = 0D;
        }
        public Ferrite(double ms, double er, double tana, double delH)
        {
            U0 = 4 * Math.PI * 1e-7;
            E0 = 8.854e-12;
            r = 2.8;
            C = 1 / Math.Sqrt(U0 * E0);
            Ms = ms; Er = er;
            losstan = tana; deltaH = delH;
        }
        #endregion
        #region [ calc ku ]
        public double Calc_ku(double Hi, double F)
        {
            ///！！！！！！！！This Function Affacts Ueff！！！！..！！！！
            double Fs = Ms * r; //固有磁饱和频率    Fs = 4πMs * r
            double F0 = Hi * r; //偏置磁场频率      F0 = H0 * r
            //工作频率          F
            double k = 0D; double u = 1D;
            double denominator = (F0 * F0) - (F * F);
            if (denominator != 0D)
            {
                u = 1 + (Fs * F0) / denominator;
                k = Fs * F / denominator;   //Fs * F  ?  F0 * F
                Ueff = (u * u - k * k) / u;
            }
            else
            {
                //谐振时的特殊情况.
                k = F / Fs;
                u = 0D;
                Ueff = 1e30;
            }
            return k / u;
        }
        public double[] Calc_ku(double[] Hi, double[] F)
        {
            double[] ku = new double[Hi.Length];
            if (Hi.Length - F.Length != 0) return ku;
            for (int i = 0; i < Hi.Length; i++)
                ku[i] = Calc_ku(Hi[i], F[i]);
            return ku;
        }
        #endregion
        #region [ calc loss]
        public double Calc_loss(double Hi, double F)
        {
            cplx A, B, C1, C2, F0;
            double Fs = Ms * r;        //固有磁饱和频率    Fs = 4πMs * r
            F0.x = Hi * r;             //偏置磁场频率      F0 = H0 * r
            F0.y = r * deltaH / 2;     //F0 = F0 + j r * deltaH/2
            A.x = 0; A.y = 1;//A = j
            B.x = 1D;// 2 * Math.PI * F;
            B.y = 0;//B = 2πf
            C1 = A * B; // C = j (2πf)
            A.x = Er; A.y = -1D * Er * losstan;                        //A = Er * (1 - losstan)            
            B = new cplx(Fs, 0) / (F0 - new cplx(F, 0));               //B = Fs/(F0-F)         
            B = new cplx(1, 0) + B;                                    //B = 1 + Fs/(F0+F)                           
            C2 = (A * B);                                                //C2= A * B = E * (u-k)                      
            C2 = C2.Pow(0.5D);
            //MessageBox.Show("C1 = " + C1.ToString() + ""n" + "C2 = " + C2.ToString() + ""n");
            return (C1 * C2).x;                       //a = j (2πf) * sqrt [E * (u+k)]
        }
        public double[] Calc_loss(double[] Hi, double[] F)
        {
            double[] loss = new double[Hi.Length];
            if (Hi.Length - F.Length != 0) return loss;
            for (int i = 0; i < Hi.Length; i++)
                loss[i] = Calc_loss(Hi[i], F[i]);
            return loss;
        }
        #endregion
        #region [ Calc Ueff ]
        public double[] Calc_Ueff(double[] Hi, double[] F)
        {
            double buffer = Ueff;       //缓存Ueff
            double[] Ue = new double[Hi.Length];
            if (Hi.Length - F.Length != 0) return Ue;
            for (int i = 0; i < Hi.Length; i++)
            {
                Calc_ku(Hi[i], F[i]);
                Ue[i] = Ueff;
            }
            Ueff = buffer;              //恢复Ueff
            return Ue;
        }
        #endregion
        #region [ Calc Lamda ]
        public double Calc_lamda(double f)
        {
            //f in MHz
            double lamda = C / (f * 1e3 * Math.Sqrt(Er * Ueff));
            return lamda;
        }
        #endregion
        #region [ Calc Y ]
        public double Calc_y(double Hi, double f, double b, double R, double a, double r)
        {
            //f in MHz, R in mm;Angles in arc
            double ku = Calc_ku(Hi, f);
            K = 2 * pi * f * 1e6 * Math.Sqrt(Er * Ueff) / C;
            double Kr = K * r * 1e-3;
            double J1 = Calc.Jn(1, Kr);
            double J0 = Calc.Jn(0, Kr);
            double J1_ = Calc.Jd(1, Kr);
            double theat = 2 * pi * (R - r) / Calc_lamda(f); //MessageBox.Show("theata = " + theat.ToString());
            double kesai_1 = 3 * b * Math.Tan(theat) / pi; //MessageBox.Show("kesai_1 = " + kesai_1.ToString());
            double kesai_2 = 3 * Math.Sin(b) * Math.Tan(theat) / pi;// MessageBox.Show("kesai_2 = " + kesai_2.ToString());
            double p = J1 / J0 + kesai_1; //MessageBox.Show("p = " + p.ToString());
            double q = J1_ / J1 - kesai_2; //MessageBox.Show("q = " + q.ToString());
            double y2 = 12 / (p * q) + 3 / (p * p); //MessageBox.Show("y2 = " + y2.ToString());
            double y = Math.Sqrt(y2); //MessageBox.Show("y = " + y.ToString());
            double M = y2 + 9 / p * p;
            double N1 = y * y * y / (Math.Sqrt(3D) * Kr);
            double N2 = Math.Sqrt(3D) * y / (p * p * Kr);
            ku_c = M / (N1 - N2);
            return y;
        }
        #endregion
        #region [ Calc Zn ]
        public double Calc_Zn(double Hi, double f, double a, double r, int n)
        {
            //f in MHz, R in mm;Angles in arc
            double ku = Calc_ku(Hi, f);
            K = 2 * pi * f * 1e6 * Math.Sqrt(Er * Ueff) / C;
            double Kr = K * r * 1e-3;
            double Jn = Calc.Jn(n, Kr);
            double Jn_ = Calc.Jd(n, Kr);
            double M;
            if (n == 0)
                M = Math.Sin(a);
            else
                M = Math.Sin(n * a) / n;
            double N = (Jn_ / Jn) - (ku * n / Kr);
            //MessageBox.Show(cplx.Exp_j(pi/3).ToString());
            return M / N;
        }
        #endregion
        #region [ Calc Zin ]
        public cplx Calc_Zin(double Hi, double f, double a, double r, double Zf, int order)
        {
            double Z0_s, Z1_s, Z2_s;
            if (order == 1)
            {
                Z0_s = Calc_Zn(Hi, f, a, r, 3) + Calc_Zn(Hi, f, a, r, -3);
                Z1_s = Calc_Zn(Hi, f, a, r, -2) + Calc_Zn(Hi, f, a, r, 4);
                Z2_s = Calc_Zn(Hi, f, a, r, 2) + Calc_Zn(Hi, f, a, r, -4);
            }
            else if (order == 2)
            {
                Z0_s = Calc_Zn(Hi, f, a, r, 6) + Calc_Zn(Hi, f, a, r, -6);
                Z1_s = Calc_Zn(Hi, f, a, r, -5) + Calc_Zn(Hi, f, a, r, 7);
                Z2_s = Calc_Zn(Hi, f, a, r, 5) + Calc_Zn(Hi, f, a, r, -7);
            }
            else
            {
                Z0_s = Calc_Zn(Hi, f, a, r, 0);
                Z1_s = Calc_Zn(Hi, f, a, r, 1);
                Z2_s = Calc_Zn(Hi, f, a, r, -1);
            }
            cplx Z0_ = new cplx(0, 3 * Zf * Z0_s / pi);
            cplx Z1_ = new cplx(0, 3 * Zf * Z1_s / pi);//Z+
            cplx Z2_ = new cplx(0, 3 * Zf * Z2_s / pi);//Z-
            double ang = 2 * pi / 3;
            cplx Z11 = (Z0_ + Z1_ + Z2_) / 3D;
            cplx Z12 = (Z0_ + Z1_ * cplx.Exp_j(ang) + Z2_ * cplx.Exp_j(0 - ang)) / 3D;
            cplx Z13 = (Z0_ + Z1_ * cplx.Exp_j(0 - ang) + Z2_ * cplx.Exp_j(ang)) / 3D;
            cplx A = Z12 * Z12 / Z13;
            return Z11 - A;
        }
        #endregion
    }
}