LPF/calcPflow.m at master · eranschweitzer/LPF · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
function P = calcPflow(ID,btype,vars,F,T,E,Sb,varargin)

P = struct();
[g,b,~,tau,tshift,~,~] = branchparts_unpack(Sb,btype);
% Gwtype = varargin_parse(varargin,'Gw','default');
% Gw = genweights(g,Gwtype);
Gw = varargin_parse(varargin,'Gw',branchweights(Sb));
switch ID
    case 0 %full
        P.f = F*vars.v.*( exp(Gw.g*F*vars.u0).*g.*tau.^(-2).*(1 + F*vars.uhat)...
                   -exp(Gw.g*T*vars.u0).*tau.^(-1).*(g.*cos(tshift) - b.*sin(tshift)).*(1 + T*vars.uhat - vars.phi)...
                   -exp(T*vars.u0).*tau.^(-1).*(g.*sin(tshift) + b.*cos(tshift)).*E*vars.theta);
        P.t = T*vars.v.*( exp(Gw.g*T*vars.u0).*g.*(1 + T*vars.uhat)...
                   -exp(Gw.g*F*vars.u0).*tau.^(-1).*( g.*cos(tshift) + b.*sin(tshift)).*(1 + F*vars.uhat - vars.phi)...
                   +exp(F*vars.u0).*tau.^(-1).*(-g.*sin(tshift) + b.*cos(tshift)).*E*vars.theta);
    case 1 %simplified
        P.f = F*vars.v.*( exp(Gw.g*F*vars.u0).*g.*(1 + F*vars.uhat)...
                   -exp(Gw.g*T*vars.u0).*(g.*cos(tshift) - b.*sin(tshift)).*(1 + T*vars.uhat - vars.phi)...
                   -exp(T*vars.u0).*(g.*sin(tshift) + b.*cos(tshift)).*E*vars.theta);
        P.t = T*vars.v.*( exp(Gw.g*T*vars.u0).*g.*(1 + T*vars.uhat)...
                   -exp(Gw.g*F*vars.u0).*( g.*cos(tshift) + b.*sin(tshift)).*(1 + F*vars.uhat - vars.phi)...
                   +exp(F*vars.u0).*(-g.*sin(tshift) + b.*cos(tshift)).*E*vars.theta);
    case 2 %no initial voltage
        P.f = F*vars.v.*( g.*(1 + F*vars.uhat)...
                   -(g.*cos(tshift) - b.*sin(tshift)).*(1 + T*vars.uhat - vars.phi)...
                   -(g.*sin(tshift) + b.*cos(tshift)).*E*vars.theta);
        P.t = T*vars.v.*( g.*(1 + T*vars.uhat)...
                   -( g.*cos(tshift) + b.*sin(tshift)).*(1 + F*vars.uhat - vars.phi)...
                   +(-g.*sin(tshift) + b.*cos(tshift)).*E*vars.theta);
    case 3 %no conductance
        P.f = F*vars.v.*(...
                   +exp(T*vars.u0).*b.*sin(tshift).*(1 + T*vars.uhat - vars.phi)...
                   -exp(T*vars.u0).*b.*cos(tshift).*E*vars.theta);
        P.t = T*vars.v.*(...
                   -exp(F*vars.u0).*b.*sin(tshift).*(1 + F*vars.uhat - vars.phi)...
                   +exp(F*vars.u0).*b.*cos(tshift).*E*vars.theta);
    case 4 %no conductance, no initial voltage
        P.f = F*vars.v.*(...
                   +b.*sin(tshift).*(1 + T*vars.uhat - vars.phi)...
                   -b.*cos(tshift).*E*vars.theta);
        P.t = T*vars.v.*(...
                   -b.*sin(tshift).*(1 + F*vars.uhat - vars.phi)...
                   +b.*cos(tshift).*E*vars.theta);
    case 5 %full nonlinear
        P.f = F*(vars.v.^2).*tau.^(-2).*g + (F*vars.v).*(T*vars.v).*...
             (cos(E*vars.theta).*(-tau.^(-1).*(g.*cos(tshift) - b.*sin(tshift)))...
             +sin(E*vars.theta).*(-tau.^(-1).*(g.*sin(tshift) + b.*cos(tshift))));
        P.t = T*(vars.v.^2).*g + (T*vars.v).*(F*vars.v).*...
             (cos(E*vars.theta).*(-tau.^(-1).*( g.*cos(tshift) + b.*sin(tshift)))...
             -sin(E*vars.theta).*(-tau.^(-1).*(-g.*sin(tshift) + b.*cos(tshift))));
    otherwise
        error('Incorrect ID')
end