classdef SuperconductingMagneticEnergyStorage < SoModel
    % Superconducting Magnetic Energy Storage (TEAM #22) Analytical Model by Ramzi Ben Ayed, 9/2011
    % S. Brisset 8/2018
    
    methods
        
        % declarative part
        
        function obj = SuperconductingMagneticEnergyStorage
            
            % create object and declare name and version
            obj = obj@SoModel('Superconducting Magnetic Energy Storage (TEAM #22) Analytical Model','by Ramzi Ben Ayed, 2011');
            
            % declare inputs of model
            obj.addInput('R1',struct('name','R1','comment','rayon bobine 1','unit','m','type','continuous'));
            obj.addInput('R2',struct('name','R2','comment','rayon bobine 2','unit','m','type','continuous'));
            obj.addInput('h1',struct('name','h1','comment','hauteur bobine 1','unit','m','type','continuous'));
            obj.addInput('h2',struct('name','h2','comment','hauteur bobine 2','unit','m','type','continuous'));
            obj.addInput('d1',struct('name','d1','comment','largeur bobine 1','unit','m','type','continuous'));
            obj.addInput('d2',struct('name','d2','comment','largeur bobine 2','unit','m','type','continuous'));
            obj.addInput('J1',struct('name','J1','comment','densité de courant bobine 1','unit','a/m2','type','continuous'));
            obj.addInput('J2',struct('name','J2','comment','densité de courant bobine 2','unit','a/m2','type','continuous'));
            obj.addInput('Bref',struct('name','Bref','comment','induction magnétique de référence','unit','T','type','constant'));
            obj.addInput('Eref',struct('name','Eref','comment','énergie stockée de référence','unit','J','type','constant'));
            obj.addInput('dipole',struct('name','dipole','comment','nombre de dipoles par bobine','unit','-','type','constant'));
            
            % declare outputs of model
            obj.addOutput('OF1',struct('name','OF1','comment','objectif sur la fuite magnétique','unit','-'));
            obj.addOutput('OF2',struct('name','OF2','comment','objectif sur l''énergie','unit','-'));
            obj.addOutput('OF',struct('name','OF','comment','OF1+OF2','unit','-'));
            obj.addOutput('E',struct('name','E','comment','énergie stockée','unit','J'));
            obj.addOutput('Bstray',struct('name','Bstray','comment','induction magnétique de fuite au carré','unit','T2'));
            obj.addOutput('Bmax1',struct('name','Bmax1','comment','induction magnétique maximale bobine 1','unit','T'));
            obj.addOutput('Bmax2',struct('name','Bmax2','comment','induction magnétique maximale bobine 2','unit','T'));
            obj.addOutput('recouvrement',struct('name','recouvrement','comment','recouvrement des bobines (<0)','unit','m'));
            obj.addOutput('quench1',struct('name','quench1','comment','condition de quench bobine 1 (<0)','unit','-'));
            obj.addOutput('quench2',struct('name','quench2','comment','condition de quench bobine 2 (<0)','unit','-'));

        end
        
        % procedural part
        
        function output = compute(obj,input)
            
            % convert input structure to variables
            fnames = fieldnames(input);
            for i=1:length(fnames)
                eval([fnames{i} '=input.' fnames{i} ';']);
            end
            
            pas_h1=dipole;
            pas_r1=dipole;
            pas_h2=dipole;
            pas_r2=dipole;
            
            tab_h11=-h1/2:h1/pas_h1:h1/2;
            tab_h22=-h2/2:h2/pas_h2:h2/2;
            tab_ray11=R1-d1/2:d1/pas_r1:R1+d1/2;
            tab_ray22=R2-d2/2:d2/pas_r2:R2+d2/2;
            
            tab_h1=(nan);tab_ray1=(nan);tab_h2=(nan);tab_ray2=(nan);
            
            for jdx=1:(length(tab_h11)-1)
                tab_h1(jdx)=tab_h11(jdx)+(tab_h11(jdx+1)-tab_h11(jdx))/2;
                tab_ray1(jdx)=tab_ray11(jdx)+(tab_ray11(jdx+1)-tab_ray11(jdx))/2;
                tab_h2(jdx)=tab_h22(jdx)+(tab_h22(jdx+1)-tab_h22(jdx))/2;
                tab_ray2(jdx)=tab_ray22(jdx)+(tab_ray22(jdx+1)-tab_ray22(jdx))/2;
            end
            
            I1=h1/pas_h1*d1/pas_r1*J1;
            I2=h2/pas_h2*d2/pas_r2*J2;
            
            %induction max dans les conducteurs primaires
            calcul_W=1;
            [B_max_cond_pr,B_M1,energy1]=obj.induction_max_cond (tab_ray1,tab_h1,tab_h1,tab_ray1,tab_h2,tab_ray2,I1,I2,calcul_W);
            B_max_cond_pr=(B_max_cond_pr);
            
            %induction max dans les conducteurs secondaires
            calcul_W=2;
            [B_max_cond_sc,B_M2,energy2]=obj.induction_max_cond (tab_ray2,tab_h2,tab_h1,tab_ray1,tab_h2,tab_ray2,I1,I2,calcul_W);
            B_max_cond_sc=(B_max_cond_sc);
            W_stockee=abs(sum(sum(energy2))+sum(sum(energy1)));
            
            %induction dans les points de l espace
            calcul_W=0;
            
            abcisses_points=0:1:10;
            ordonees_points=0:1:10;
            vect_B_point1=(nan);
            for jdx=1:length(abcisses_points);
                [B_max,B_M1,WW1]=obj.induction_max_cond (abcisses_points(jdx),ordonees_points(11),tab_h1,tab_ray1,tab_h2,tab_ray2,I1,I2,calcul_W);
                vect_B_point1(jdx)=B_M1;
                B_stray1=sum(vect_B_point1)/length(vect_B_point1);
            end
            
            vect_B_point2=(nan);
            for jdx=1:length(ordonees_points);
                [B_max,B_M2,WW2]=obj.induction_max_cond (abcisses_points(11),ordonees_points(jdx),tab_h1,tab_ray1,tab_h2,tab_ray2,I1,I2,calcul_W);
                vect_B_point2(jdx)=B_M2;
                B_stray2=sum(vect_B_point2)/length(vect_B_point2);
            end
            B_stray=B_stray1+B_stray2;
            B_st_carre=B_stray^2;
            
            % outputs
            output.E = W_stockee;
            output.Bstray = B_st_carre;
            output.OF1 = B_st_carre/Bref^2;
            output.OF2 = abs(W_stockee-Eref)/Eref;
            output.OF = output.OF1+output.OF2;
            output.recouvrement = -((R2-d2/2)-(R1+d1/2));
            output.Bmax1 = B_max_cond_pr;
            output.Bmax2 = B_max_cond_sc;
            output.quench1 = 1e-6*abs(J1)+6.4*abs(output.Bmax1)-54;
            output.quench2 = 1e-6*abs(J2)+6.4*abs(output.Bmax2)-54;
            
        end
        
        % fonctions internes
        
        function [Bmax_cond,B_M,vect_W]= induction_max_cond (obj,ro_M,z_M,z_sp1,r_sp1,z_sp2,r_sp2,dI1,dI2,drapeau)
            
            vect_Bz=(nan);
            vect_W=(nan);
            tab_h1=z_sp1;
            tab_ray1=r_sp1;
            tab_h2=z_sp2;
            tab_ray2=r_sp2;
            
            for j1=1:length(ro_M)
                for j2=1:length(z_M)
                    
                    som_A1=0;
                    som_Br1=0;
                    som_Bz1=0;
                    
                    for idx=1:length(tab_h1)
                        for idx1=1:length(tab_ray1)
                            [A_phi1,Bro1,Bz1]=obj.problem22 (ro_M(j1),z_M(j2),tab_h1(idx),tab_ray1(idx1),dI1);
                            som_A1=som_A1+(A_phi1);
                            som_Br1=som_Br1+Bro1;
                            som_Bz1=som_Bz1+Bz1;
                        end
                    end
                    
                    som_A2=0;
                    som_Br2=0;
                    som_Bz2=0;
                    
                    for idx=1:length(tab_h2)
                        for idx1=1:length(tab_ray2)
                            [A_phi2,Bro2,Bz2]=obj.problem22 (ro_M(j1),z_M(j2),tab_h2(idx),tab_ray2(idx1),dI2);
                            som_A2=som_A2+(A_phi2);
                            som_Br2=som_Br2+Bro2;
                            som_Bz2=som_Bz2+Bz2;
                        end
                    end
                    
                    A_phi_M=(som_A2)+(som_A1);
                    Bz_M=som_Bz2+som_Bz1;
                    Br_M=som_Br2+som_Br1;
                    B_M=sqrt(Br_M^2+Bz_M^2);
                    vect_Bz(j1,j2)=B_M;

                    if drapeau==1
                        vect_W1(j1,j2)=pi*A_phi_M*dI1*ro_M(j1);
                        vect_W=vect_W1;
                    elseif drapeau==2
                        vect_W2(j1,j2)=pi*A_phi_M*dI2*ro_M(j1);
                        vect_W=vect_W2;
                    end
                end
            end
            
            max_vect=max(vect_Bz);
            min_vect=min(vect_Bz);
            induction_max=max(max_vect);
            induction_min=min(min_vect);
            Bmax_cond=max(abs(induction_max),abs(induction_min));
        end
        
        function [A_phi,Bro,Bz]=problem22(obj,ro,z,d,a,I)
            
            mu0=4*pi*1e-7;
            
            if (ro==a) && ((z-d)==0)% induction dans le conducteur de la spire
                A_phi=0;
                B_spire=0;
                Bz=B_spire;
                Bro=0;
            elseif (ro==0)% indcution sur l axe de la spire
                Bro=0;
                Bz=mu0*I*a^2/(2*((z-d)^2+a^2)^1.5);
                A_phi=0;
            elseif (ro>=10)||(z>=10) % loin de la spire
                A_phi=mu0*I/4*a^2*ro/(ro^2+z^2)^2.5;
                Bz=(mu0*I*a^2/4)*(2/(ro^2+(z-d)^2)^1.5-3*ro^2/(ro^2+(z-d)^2)^2.5);
                Bro=(mu0*I*a^2/4)*3*ro*(z-d)/(ro^2+(z-d)^2)^2.5;
                
            else %dans la spire
                
                k=2*sqrt(a*ro)/sqrt((a+ro)^2+(z-d)^2);
                [K_k,E_k] = ellipke(k);
                
                Bro=mu0*I*k*(z-d)/(4*pi*ro*sqrt(a*ro))*(K_k+(a^2+ro^2+(z-d)^2)/((a-ro)^2+(z-d)^2)*E_k);
                Bz=mu0*I*k/(4*pi*sqrt(a*ro))*(K_k+(a^2-ro^2-(z-d)^2)/((a-ro)^2+(z-d)^2)*E_k);
                A_phi=mu0*I/(k*pi)*sqrt(a/ro)*((1-k^2/2)*K_k-E_k);
                
            end
        end
        
    end
end