/* Соглашения о нумерации узлов.
nd - количество диодов.
1..nd+1 - последовательно пронумерованные узлы диодной цепи.
1 - первый узел диодной цепи, является общим проводом
в случае чётного коэффициента; подключается к источнику в случае нечётного.
nd+1 - последний, выходной узел диодной цепи, к нему подключается нагрузка.
Номер узла, соответствующего общему проводу, должен обязательно
быть равен 1 (чётный коэффициент), либо иметь наибольший номер
среди всех узлов.
Первый вывод источника питания подключается ко входу умножителя:
к узлу 1, если умножитель имеет нечётный коэффициент или
к некоторому конденсатору, если коэффициент чётный.
Второй вывод источника питания подключается всегда к общему проводу
(так можно определить номер узла, соответствующего общему проводу).
Если общий провод имеет номер 1, то максимальный номер должен быть у узла-
входа умножителя (к которому подключён источник).
В результате принятых соглашений, питающий источник оказывается всегда
включённым между узлами 1 и узлом с наибольшим номером.
*/

/* Соглашения о структуре матрицы соединений circuit.
Если схема умножителя содержит nc конденсаторов, матрица
будет иметь размер (nc+1)*2.
Первые строки 1..nc описывают подключение соответствующих
конденсаторов: circuit[i, 1] - номер узла, к которому подключён
первый вывод конденсатора Ci; circuit[i, 2] - номер узла, к которому
подключён второй вывод.
Последняя строка описывает подключение источника питания.
Диоды и нагрузка не указываются в матрице, способ их подключения 
является предопределённым. */

/* Присвоить элементам строки i матрицы A значения из списка (матрицы) r */

set_row(A, i, r):=block(
    [],
    if not listp(r) then r:list_matrix_entries(r),
    for j:1 thru length(r) do A[i, j]:r[j],
    A)$

/* Присвоить элементам столбца j матрицы A значения из списка (матрицы) c */

set_col(A, j, c):=block(
    [cl],
    cl:if not listp(c) then list_matrix_entries(c)
    else c,
    for i:1 thru length(cl) do A[i, j]:cl[i],
    A)$

/* Функция возвращает матрицу соединений для умножителя на nd типа A. */

create_circuitA(nd):=block(
    [n, m, circuit],
    n:nd+2,
    m:nd+1,
    circuit:zeromatrix(m, 2),
    for i:1 thru nd do (
        circuit[i, 1]:i+1,
        circuit[i, 2]:i-1
    ),
    circuit[1, 2]:n,
    circuit[m, 1]:if mod(nd, 2)=0 then n else 1,
    circuit[m, 2]:if mod(nd, 2)=0 then 1 else n,
    circuit)$

/* Функция возвращает матрицу соединений для умножителя на nd типа B. */

create_circuitB(nd):=block(
    [n, m, circuit],
    if nd<=2 then return(create_circuitA(nd)),
    n:nd+2,
    m:nd+1,
    circuit:zeromatrix(m, 2),
    for i:1 thru m do (
        circuit[i, 1]:i+1,
        circuit[i, 2]:i-3
    ),
    circuit[1, 2]:n,
    circuit[2, 2]:1,
    circuit[3, 2]:n,
    circuit[m, 2]:1,
    if mod(nd, 2)=1 then (
        circuit[m, 1]:1,
        circuit[m, 2]:n
    ),
    circuit)$

/* Функция возвращает матрицу соединений для умножителя на nd типа Z. */

create_circuitZ(nd):=block(
    [n, m, circuit],
    n:nd+2,
    m:nd+1,
    circuit:zeromatrix(m, 2),
    for i:1 thru m do (
        circuit[i, 1]:i+1,
        circuit[i, 2]:if mod(i, 2)=1 then n else 1
    ),
    circuit[m, 2]:1,
    if mod(nd, 2)=1 then (
        circuit[m, 1]:1,
        circuit[m, 2]:n
    ),
    circuit)$

/* По матрице соединений сформировать узловую матрицу,
считая диоды разомкнутыми ветвями (т.е. игнорируя их).
Автоматически добавляется ветвь - нагрузка. */

get_node_matrix(circuit, nd):=block(
    [nn, nbr, gnd_node, A], 
    
    /* TODO: можно добавить проверку валидности circuit */
    
    /* Ищем узел с наибольшим номером, считаем этот номер равным
    общему количеству узлов. */
    nn:lmax(list_matrix_entries(circuit)),
    /* Количество ветвей равно количеству строк в circuit 
     + 1 ветвь - нагрузка */  
    nbr:length(circuit)+1,
    A:zeromatrix(nn, nbr),
    for i:1 thru nbr-1 do (
        A[circuit[i, 1], i]:1,
        A[circuit[i, 2], i]:-1),
    gnd_node:circuit[nbr-1, 2],
    A[nd+1, nbr]:1,
    A[gnd_node, nbr]:-1,
    A)$

/* Возвращает узловую матрицу, полученную
из узловой матрицы A исключением висящих ветвей. */

exclude_orphan_branches(A):=block(
    [A:copymatrix(A), sz, r, k, j, n:1],
    sz:matrix_size(A),
    /* Выполняем в цикле, т.к. исключение одних висящих ветвей
    может сделать висящими другие ветви. */
    while n#0 do (
        n:0,
        for i:1 thru sz[1] do (
            r:abs(row(A, i)),
            if sum(r[1, k], k, 1, sz[2])=1 then (
                j:locate_matrix_entry(r, 1, 1, 1, sz[2], identity, 'max)[2],
                set_col(A, j, zeromatrix(sz[1], 1)),
                n:n+1))),
    A)$

/* Возвращает узловую матрицу дерева, построенного для
схемы, которая задана узловой матрицей A. */ 

get_tree_matrix(A):=block(
    [T, w, sz, c, z, ntr, Test, s],
    T:zerofor(A),
    sz:matrix_size(A),
    w:zeromatrix(sz[1], 1),
    ntr:{},
    for j:1 thru sz[2] do (
        c:abs(col(A, j)),
        z:w*c,
        s:sum(z[i, 1], i, 1, length(z)),
        if s<2 then (
            w:w+c-z,
            set_col(T, j, col(A, j))
        ) elseif s=2 then (
            Test:copymatrix(T),
            set_col(Test, j, col(A, j)),
            if zeromatrixp(exclude_orphan_branches(Test)) then (
                set_col(T, j, col(A, j))
            )
            else ntr:adjoin(j, ntr)
        )
    ),
    [T, ntr])$

/* Возвращает строку матрицы контуров, соответствующую одному контуру,
найденному по узловой матрице A.
Возвращаемый вектор-строка имеет тип список. */

get_closed_loop(A):=block(
    [A:copymatrix(A), n, m, sz, result, nl, ie, i, j, r],
    [n, m]:matrix_size(A),
    nl:makelist(0, m),
    result:copylist(nl),
    A:exclude_orphan_branches(A),
    [i, j]:locate_matrix_entry(A, 1, 1, n, m, lambda([x], x#0), 'bool),
    if i=false then return(nl),
    result[ j]:A[i, j],
    A[i, j]:0,
    [ie, j]:locate_matrix_entry(A, 1, j, n, j, lambda([x], x#0), 'bool),
    if ie=false then return(nl),
    A[ie, j]:0,
    r:while i#ie do (
        [i, j]:locate_matrix_entry(A, i, 1, i, m, lambda([x], x#0), 'bool),
        if j=false then return(false),
        result[j]:-A[i, j],
        A[i, j]:0,
        [i, j]:locate_matrix_entry(A, 1, j, n, j, lambda([x], x#0), 'bool),
        if j=false then return(false),
        A[i, j]:0
    ),
    if r=false then return(nl),
    result)$

/* Возвращает матрицу контуров для схемы, заданной узловой матрицей A. */

get_independent_loops(A):=block(
    [T, ntr, lst, result, r, tmp],
    [T, ntr]:get_tree_matrix(A),
    lst:listify(ntr),
    result:matrix(),
    for j:1 thru length(lst) do (
        tmp:copymatrix(T),
        set_col(tmp, lst[j], col(A, lst[j])),
        r:get_closed_loop(tmp),
        if not zeromatrixp(r) then result:addrow(result, r)
    ),
    result)$

/* Возвращает узловую матрицу, получаемую из A при условии
замыкания нечётных (k=1) или чётных (k=2) диодов. */

switch_diodes(A, nd, k):=block(
    [result, n, m, n2, nj, r1, r2, di, i0],
    if k#1 and k#2 then error("Недопустимое значение аргумента k."),
    [n, m]:matrix_size(A),
    nj:if mod(nd, 2)=1 then (nd+3-2*k)/2 else nd/2,
    n2:n-nj,
    result:zeromatrix(n2, m),
    di:if k=1 then 0 else 1,
    if di>0 then set_row(result, 1, row(A, 1)),
    for i:1 thru nj do (
        r1:row(A, 2*i-1+di),
        r2:row(A, 2*i+di),
        if not zeromatrixp(r1*r2) then error("Недопустимая узловая матрица: конденсатор включён между соседними узлами диодной цепи."),
        set_row(result, i+di, r1+r2)
    ),
    i0:2*nj+di+1,
    for i:i0 thru n do set_row(result, i-nj, row(A, i)),
    result)$

/* Возвращает уравнения переключения - уравнения цепи (в виде списка уравнений),
составленные по узловой матрице A при условии замыкания нечётных (k=1) или
чётных диодов (k=2).
nd - количество диодов в умножителе.
x1, x2 - список переменных в начальный и конечный моменты времени.
xc - список с ёмкостями конденсаторов. */

get_switch_equations(A, nd, k, x1, x2, xc):=block(
    [As, An, Y, X, eq, n, m, Q],
    As:switch_diodes(A, nd, k),
    
    /* Узловые уравнения */
    An:submatrix(1, length(As), As),
    [n, m]:matrix_size(An),
    X:if m=0 then [] else append(makelist(xc[i]*(x2[i]-x1[i]), i, 1, m-2), [0, 0]),
    Y:An.X,
    if not matrixp(Y) then Y:matrix([Y]),
    eq:makelist(Y[i, 1]=0, i, 1, n),
    
    /* Контурные уравнения */
    Q:get_independent_loops(As),
    [n, m]:matrix_size(Q),
    if n=0 then return(eq),
    X:makelist(x2[i], i, 1, m),
    Y:Q.X,
    if not matrixp(Y) then Y:matrix([Y]),
    n:length(Q),
    eq:append(eq, makelist(Y[i, 1]=0, i, 1, n)),
    eq)$

/* Функция составляет систему уравнений для заданного умножителя и решает её.
Возвращает список из трёх элементов: выходное напряжение, размах пульсаций и
список из всех найденных решений системы уравнений (напряжения на конденсаторах,
мгновенное напряжение на источнике и на выходе в характерных точках). */

solve_multiplier(circuit, nd, xc, a, F, IL):=block(
    [u, ux, u1, uy, u2,
        A, T, nbr, n, m, out_node, gnd_node, power_node,
        M1, M2, X1, X2, Y, eq, g, k, xlst, r, result, uout, hout],
    if not listp(xc) or length(xc)#length(circuit)-1 then
        error("Аргумент xc должен быть списком, количество элементов которого равно количеству конденсаторов в схеме."),
    A:get_node_matrix(circuit, nd),
    [n, m]:matrix_size(A),
    [T, nbr]:get_tree_matrix(A),
    out_node:nd+1,
    gnd_node:circuit[length(circuit), 2],
    power_node:circuit[length(circuit), 1],
    eq:[],
    
    /* Уравнения разряда при закрытых диодах, 2 однотипных системы. */
    M1:submatrix(gnd_node, power_node, A),
    X1:makelist(0, m),
    X2:makelist(0, m),
    for i:1 thru m-2 do (
        X1[i]:xc[i]*(ux[i]-u[i]),
        X2[i]:xc[i]*(uy[i]-u1[i])
    ),
    X1[m]:IL/(2*F),
    X2[m]:IL/(2*F),
    Y:M1.X1,
    if not matrixp(Y) then Y:matrix([Y]),
    /* Получаем уравнений на 2 меньше, чем узлов */
    eq:append(eq, makelist(Y[i, 1]=0, i, 1, n-2)),
    Y:M1.X2,
    if not matrixp(Y) then Y:matrix([Y]),
    eq:append(eq, makelist(Y[i, 1]=0, i, 1, n-2)),
    /* Если есть замкнутые контуры, добавляем уравнения для них */
    M2:get_independent_loops(A),
    g:length(M2),
    if g>0 then (
        X1:makelist(ux[i], i, 1, m),
        X2:makelist(uy[i], i, 1, m),
        Y:M2.X1, if not matrixp(Y) then Y:matrix([Y]),
        eq:append(eq, makelist(Y[i, 1]=0, i, 1, g)),
        Y:M2.X2, if not matrixp(Y) then Y:matrix([Y]),
        eq:append(eq, makelist(Y[i, 1]=0, i, 1, g))
    ),
    
    /* Уравнения первого перезаряда: ux->u1 */
    k:if mod(nd, 2)=0 then 1 else 2,
    X1:makelist(ux[i], i, 1, m),
    X2:makelist(u1[i], i, 1, m),
    eq:append(eq, get_switch_equations(A, nd, k, X1, X2, xc)),
  
    /* Уравнения второго перезаряда: uy->u2 */
    k:if mod(nd, 2)=0 then 2 else 1,
    X1:makelist(uy[i], i, 1, m),
    X2:makelist(u2[i], i, 1, m),
    eq:append(eq, get_switch_equations(A, nd, k, X1, X2, xc)),
    
    /* Уравнения стационарного режима: u[i]=u2[i] */
    eq:append(eq, makelist(u[i]=u2[i], i, 1, m)),
    
    /* Уравнения питающего умножитель источника переменного напряжения. */
    eq:append(eq, [u[m-1]=a, ux[m-1]=-a, u1[m-1]=-a, uy[m-1]=a]),
    
    /* Решаем полученную систему относительно следующих переменных: */
    xlst:append(
        makelist(u[i], i, 1, m),
        makelist(ux[i], i, 1, m),
        makelist(u1[i], i, 1, m),
        makelist(uy[i], i, 1, m),
        makelist(u2[i], i, 1, m)),
    r:solve(eq, xlst),
    if length(r)=0 then error("Ошибка: не найдено решение системы."),
    /* Преобразуем результат решения в список значений переменных
    в порядке их перечисления в xlst, значения u2[i] не включаем в список */
    result:makelist(args(r[1][i])[2], i, 1, 4*m),
    uout:result[m], hout:uout-result[4*m],
    [expand(uout), expand(hout), expand(result)]
)$

/* Тестирование функции solve_multiplier на заведомо известных решениях */
test_solve_multiplier(iterations):=block(
    [tp:["A", "B", "Z"], type, create_circuit, err:1/100,
        d, f, fh, fa, fha, fb, fhb, fz, fhz, uout, hout, r, th_uout, th_hout, cir, result:"OK"],
    d:IL/(2*F*C),
    fa(m):=m*a-d*(
        if mod(m, 2)=0 then m^3/6+m^2/8+m/12
        else m^3/6+m^2/8-m/6-1/8),
    fha(m):=d*if mod(m, 2)=0 then (m^2/4+m/2) else (m^2/4+m+3/4),
    fb(m):=m*a-d*(m^3/24+5*m^2/32+
        if mod(m, 4)=0 then 5*m/24
        elseif mod(m, 4)=1 then m/48-7/32
        elseif mod(m, 4)=2 then m/3+3/8
        else /*mod(m, 4)=3*/ 19*m/48+9/32),
    fhb(m):=d*(m^2/16+
        if mod(m, 4)=0 then m/4
        elseif mod(m, 4)=1 then 5*m/8+21/16
        elseif mod(m, 4)=2 then m/2+3/4
        else 3*m/8+5/16),
    fz(m):=m*a-d*(2*m-2),
    fhz(m):=2*d,
    f(type, m):=if type="A" then fa(m)
        elseif type="B" then fb(m)
        elseif type="Z" then fz(m)
        else error(type, "is an unknown type of circuits."),
    fh(type, m):=if type="A" then fha(m)
        elseif type="B" then fhb(m)
        elseif type="Z" then fhz(m)
        else error(type, "is an unknown type of circuits."),
    create_circuit(type, m):=if type="A" then create_circuitA(m)
        elseif type="B" then create_circuitB(m)
        elseif type="Z" then create_circuitZ(m)
        else error(type, "is an unknown type of circuits."),
    
    prederror:true,
     
    display("Type A, x12, C1=4C"),
    th_uout:12*a-253*d,
    th_hout:42*d,
    display(th_uout, th_hout),
    xc:makelist(C, 12),
    xc[1]:4*C,
    [uout, hout, r]:solve_multiplier(create_circuitA(12), 12, xc, a, F, IL),
    if not equal(th_uout, uout) or not equal(th_hout, hout) then error("Обнаружено расхождение."),
    display(result),
    
    display("Type Design1, x4"),
    th_uout:4*a-17/3*d,
    th_hout:2*d,
    display(th_uout, th_hout),
    xc:makelist(C, 6),
    cir:matrix([2, 6], [3, 1], [4, 2], [5, 3], [4, 6], [5, 1], [6, 1]), 
    [uout, hout, r]:solve_multiplier(cir, 4, xc, a, F, IL),
    if not equal(th_uout, uout) or not equal(th_hout, hout) then error("Обнаружено расхождение."),
    display(result),
    
    for nd:1 thru iterations do (
        for j:1 thru length(tp) do (
            type:tp[j],
            th_uout:f(type, nd),
            th_hout:fh(type, nd),
            xc:makelist(C, nd),
            display(type, nd, th_uout, th_hout),
            [uout, hout, r]:solve_multiplier(create_circuit(type, nd), nd, xc, a, F, IL),
            display(uout, hout),
            if not equal(th_uout, uout) or not equal(th_hout, hout)
            then error("Обнаружено расхождение."),
            display(result)
        )
    ),
    true)$