#!/usr/bin/env python
# -*- coding: utf-8 -*-
u'''
Created on 2. 5. 2014
Komplet celý výpočet dle normy ČSN EN 13445-3
@author: Stanislav Schütz
'''
def Soucinitele(g_0,g_1,h,l_0,mi):
    A = g_1/g_0 - 1
    C = 48*(1-mi*mi)*(h/l_0)**4
    C_1 = 1.0/3 + A/12.0
    C_2 = 5.0/42 + 17*A/336
    C_3 = 1.0/210 + A/360
    C_4 = 11.0/360 + 59*A/5040 + (1+3*A)/C
    C_5 = 1.0/90 + 5*A/1008 - (1 + A)**3/C
    C_6 = 1.0/120 + 17*A/5040 + 1/C
    C_7 = 215.0/2772 + 51*A/1232 + (120+225*A+150*A*A+35*A**3)/14/C
    C_8 = 31.0/6930 + 128*A/45045 + (66+165*A+132*A*A+35*A**3)/77/C
    C_9 = 533.0/30240 + 653*A/73920 + (42+198*A+117*A*A+25*A**3)/84/C
    C_10 = 29.0/3780 + 3*A/704 - (42+198*A+243*A*A+91*A**3)/84/C
    C_11 = 31.0/6048 + 1763*A/665280 + (42+72*A+45*A*A+10*A**3)/84/C
    C_12 = 1.0/2925 + 71*A/300300 + (88+198*A+156*A*A+42*A**3)/385/C
    C_13 = 761.0/831600 + 937*A/1663200 + (2+12*A+11*A*A+3*A**3)/70/C
    C_14 = 197.0/415800 + 103*A/332640 - (2+12*A+17*A*A+7*A**3)/70/C
    C_15 = 233.0/831600 + 97*A/554400 + (6+18*A+15*A*A+4*A**3)/210/C
    C_16 = C_1*C_7*C_12 + C_2*C_8*C_3 + C_3*C_8*C_2 - (C_3*C_3*C_7 + C_8*C_8*C_1 + C_2*C_2*C_12) # Pozor!! Dva stejné členy!!
    C_17 = (C_4*C_7*C_12 + C_2*C_8*C_13 + C_3*C_8*C_9 - (C_13*C_7*C_3 + C_8*C_8*C_4 + C_12*C_2*C_9))/C_16
    C_18 = (C_5*C_7*C_12 + C_2*C_8*C_14 + C_3*C_8*C_10 - (C_3*C_7*C_14 + C_8*C_8*C_5 + C_12*C_2*C_10))/C_16
    C_19 = (C_6*C_7*C_12 + C_2*C_8*C_15 + C_3*C_8*C_11 - (C_3*C_7*C_15 + C_8*C_8*C_6 + C_12*C_2*C_11))/C_16
    C_20 = (C_1*C_9*C_12 + C_4*C_8*C_3 + C_3*C_13*C_2 - (C_3*C_3*C_9 + C_13*C_8*C_1 + C_12*C_4*C_2))/C_16
    C_21 = (C_1*C_10*C_12 + C_5*C_8*C_3 + C_3*C_14*C_2 - (C_3*C_3*C_10 + C_14*C_8*C_1 + C_12*C_5*C_2))/C_16
    C_22 = (C_1*C_11*C_12 + C_6*C_8*C_3 + C_3*C_15*C_2 - (C_3*C_3*C_11 + C_15*C_8*C_1 + C_12*C_6*C_2))/C_16
    C_23 = (C_1*C_7*C_13 + C_2*C_9*C_3 + C_4*C_8*C_2 - (C_3*C_7*C_4 + C_8*C_9*C_1 + C_2*C_2*C_13))/C_16
    C_24 = (C_1*C_7*C_14 + C_2*C_10*C_3 + C_5*C_8*C_2 - (C_3*C_7*C_5 + C_8*C_10*C_1 + C_2*C_2*C_14))/C_16
    C_25 = (C_1*C_7*C_15 + C_2*C_11*C_3 + C_6*C_8*C_2 - (C_3*C_7*C_6 + C_8*C_11*C_1 + C_2*C_2*C_15))/C_16
    C_26 = -(C/4)**0.25
    C_27 = C_20 - C_17 - 5.0/12 + C_17*C_26
    C_28 = C_22 - C_19 - 1.0/12 + C_19*C_26
    C_29 = -(C/4)**0.5
    C_30 = -(C/4)**0.75
    C_31 = 3*A/2 - C_17*C_30
    C_32 = 0.5 - C_19*C_30
    C_33 = C_26*C_32/2 + C_28*C_31*C_29 - (C_30*C_28/2 + C_32*C_27*C_29)
    C_34 = 1.0/12 + C_18 - C_21 - C_18*C_26
    C_35 = C_18*C_30
    C_36 = (C_28*C_35*C_29 - C_32*C_34*C_29)/C_33
    C_37 = (C_26*C_35/2 + C_34*C_31*C_29 - C_34*C_30/2 - C_35*C_27*C_29)/C_33
    
    E_1 = C_17*C_36 + C_18 + C_19*C_37
    E_2 = C_20*C_36 + C_21 + C_22*C_37
    E_3 = C_23*C_36 + C_24 + C_25*C_37
    E_4 = (3 + C_37 + 3*C_36)/12 - (2*E_3 + 15*E_2 + 10*E_1)/10
    E_5 = E_1*(3+A)/6 + E_2*(21+11*A)/84 + E_3*(3+2*A)/210
    E_6 = E_5 - C_36*(7.0/120+A/36+3*A/C) - 1.0/40 - A/72 - C_37*(1.0/60+A/120+1/C)
    
    souc = {'A':A,          # slovník určený pro výstup
            'C':C,
            'C_18':C_18,
            'C_21':C_21,
            'C_24':C_24,
            'C_36':C_36,
            'E_4':E_4,
            'E_6':E_6}
    return souc


def Srouby(typ,w,n,d,g_1,B,C,G,P,f_BA,f_B,y,m):
    import math
        
    b_0 = w/2
    if b_0 <= 6.3: # podmínka, která spočte účinnou šířku těsnění
        b = b_0    
    else:
        b = 2.52 * math.sqrt(b_0)
        G = G + w/2 - 2*b # ze středního průměru těsnění vypočte působiště reakční síly v těsnění
    
    H = math.pi / 4 * G*G*P # H je osová síla od tlaku
    
    H_G = 2* math.pi * G *b*m*P # tlaková síla pro dasažení těsnosti
    H_D = B*B*P*math.pi/4 # osová síla od tlaku přenášená přes plášt na přírubu
    H_T = H - H_D # osová síla od tlaku na čelní plochu příruby
    
    if typ == 1:
        h_D = (C-B-g_1)/2 # vzdálenost mezi šroubem a působištěm H_D v závislosti na typu
    else:
        h_D = (C-B)/2
    h_G = (C-G)/2 #vzdálenost mezi šroubem a působištěm reakční síly těsnění H_G
    h_T = (2*C-B-G)/4 # vzdálenost mezi šroubem a působištěm H_T
    
    # síly ve šroubech:
    W_a = G*b*y* math.pi # min požadovaná síla při MONTÁŽI    
    W_op = H + H_G # min požadovaná síla při PROVOZU
    
    # průřezy šroubů
    A_Bmin = max(W_a/f_BA, W_op/f_B) # požadovaný průřez šroubů
    A_B = n*d*d*0.25*math.pi # skutečný průřez šroubů
    
    # výsledná síla a momenty
    W = 0.5*(A_Bmin+A_B)*f_BA # konstrukční síla ve šroubech při montáži, kde beru P=0, tedy Amin = 0
    #print W
    M_A = W*h_G # celkový montážní moment
    M_op = H_D*h_D + H_T*h_T + H_G*h_G # provozní moment
    
    zatizeni = {'M_A':M_A,
                'M_op':M_op,
                'A_B':A_B,
                'A_Bmin':A_Bmin,
                'W':W}
    return zatizeni


def Mat(mat,t,Tep):
    # P235GH
    mat1 = [(235,227,214,198,182,167,153,142,133), # t<=16
            (225,218,205,190,174,160,147,136,128), # t<=40
            (215,208,196,181,167,153,140,130,122), # t<=60
            (200,193,182,169,155,142,130,121,114), # t<=100
            (185,179,168,156,143,131,121,112,105), # t<=150
            360.0]
    # P355GH
    mat2 = [(355,343,323,299,275,252,232,214,202),
            (345,334,314,291,267,245,225,208,196),
            (335,324,305,282,259,238,219,202,190),
            (315,305,287,265,244,224,206,190,179),
            (295,285,268,249,228,209,192,178,167),
            510.0]
    # 20MnMoNi4-5
    mat3 = [(470,460,448,439,432,424,415,402,384),
            (470,460,448,439,432,424,415,402,384),
            (460,450,438,430,423,415,406,394,375),
            (450,441,429,420,413,406,398,385,367),
            (440,431,419,411,404,397,389,377,359),
            (590.0,570.0)]
    # 12CrMo9-10
    mat4 = [(355,341,323,311,303,298,295,292,287),
            (355,341,323,311,303,298,295,292,287),
            (355,341,323,311,303,298,295,292,287),
            (355,341,323,311,303,298,295,292,287),
            (355,341,323,311,303,298,295,292,287),
            540]
    
    # Určí, do kterého intervalu tlouštka spadá - přiřazené číslo odpovídá indexu řádku
    if t <= 16:
        ind = 0
    elif t<= 40:
        ind = 1
    elif t<= 60:
        ind = 2
    elif t<= 100:
        ind = 3
    elif t<= 150:
        ind = 4
    
    # Výběr materiálu provedeme formou čísla    
    if mat == 1:
        mat = mat1
    elif mat == 2:
        mat = mat2
    elif mat == 3:
        mat = mat3
    elif mat == 4:
        mat = mat4
    
    hodnoty = mat[ind] # Sem vezmu řádek ze zvoleného materiálu a odpovídající tloušťky
        
    index = int(Tep/50)
    
    # Lineární interpolace dle Teploty, navrácení hodnoty Re
    if index == 0 and Tep >= 20:
        d = hodnoty[index]
        a = float(Tep - 20)/30
        b = float(hodnoty[index+1]-hodnoty[index])
        c = a*b
        Re = d + c
    elif Tep < 20:
        Re = float(hodnoty[0]) # pro teploty menší než 20 bereme konzervativní hodnotu pro T = 20°C
    elif Tep <= 400:
        d = hodnoty[index]
        a = float(Tep - index*50)/50
        b = float(hodnoty[index+1]-hodnoty[index])
        c = a*b
        Re = d + c
    
    # Určení hodnoty meze pevnosti
    if mat == mat3: # Vyjímka pro mat3 - vliv tloušťky na Rm 
        rozmezi = mat[-1]
        if t <= 60:
            Rm = rozmezi[0]
            Rm = float(Rm)
        else:
            Rm = rozmezi[1]
            Rm = float(Rm)
    else:
        Rm = mat[-1]    
        
    f = min(Re/1.5,Rm/2.4) # Dovolené namáhání dle normy
    return {'f':f,
            'Re':Re,
            'Rm':Rm}
   
    
def volna(e,betaY,M):
    u'''
    Created on 13. 3. 2014
    11.5.4.1.1 - Volná metoda
    Funkce vrátí hodnoty napětí podle volné metody
    e - nejmenší tlouštka příruby
    M - moment pro dané zatížení
    betaY - koeficient
    @author: Stanislav Schütz
    '''
    # Výpočet jednotlivých napětí:
    sigmaTangenc = betaY*M/e/e
    sigmaR = 0
    sigmaH = 0

    sigma = {'sigmaR':sigmaR,
             'sigmaTangenc':sigmaTangenc,
             'sigmaH':sigmaH}
    return sigma


def integralne(g_0,g_1,e,h,K,l_0,betaT,betaU,betaY,mi,M):
    u'''
    Created on 16. 3. 2014
    11.5.4.1.1 - Integrální metoda
    Vstupy pro soucinitele viz soucinitele.py
    e - min tlouštka příruby
    K, betaT,U,Y - součinitele pro výpočet
    M - zatěžující moment, nutno ještě pořešit
    @author: Stanislav Schütz
    '''
    souc = Soucinitele(g_0,g_1,h,l_0,mi) # vloženi součinitelů
    
    betaF = -souc['E_6']/(souc['C']/3/(1-mi*mi))**0.25/(1+souc['A'])**3*souc['C']
    betaV = souc['E_4']/(3*(1-mi*mi)/souc['C'])**0.25/(1+souc['A'])**3
    fix = souc['C_36']/(1+souc['A'])
    fi = max(fix,1)
    lamb = (e*betaF+l_0)/betaT/l_0 + e*e*e*betaV/betaU/l_0/g_0/g_0
    
    sigmaH = fi*M/lamb/g_1/g_1
    sigmaR = (1.333*e*betaF+l_0)*M/lamb/e/e/l_0
    sigmaTangenc = betaY*M/e/e - sigmaR*(K*K+1)/(K*K-1)
    sigma = {'sigmaR':sigmaR,
             'sigmaTangenc':sigmaTangenc,
             'sigmaH':sigmaH}
    return sigma


def skrkem(g_0,g_1,e,h,K,l_0,betaT,betaU,betaY,mi,M):
    u'''
    Created on 17. 3. 2014
    11.5.4.1.1 - Volná metoda přírub s krkem
    Vstupy pro soucinitele viz soucinitele.py
    e - min tlouštka příruby
    K, betaT,U,Y - součinitele pro výpočet
    M - zatěžující moment, nutno ještě pořešit
    @author: Stanislav Schütz
    '''
    souc = Soucinitele(g_0,g_1,h,l_0,mi) # vloženi součinitelů
    
    betaFL = -(souc['C_18']*(3+souc['A'])/6 + souc['C_21']*(21+11*souc['A'])/84 + 
              souc['C_24']*(3+2*souc['A'])/210 - 
              (9+5*souc['A'])/360)/(souc['C']/3/(1-mi*mi))**0.25/(1+souc['A'])**3*souc['C']
    betaVL = (0.25 - souc['C_24']/5 - 1.5*souc['C_21'] - 
              souc['C_18'])/(3*(1-mi*mi)/souc['C'])**0.25/(1+souc['A'])**3
    lamb = (e*betaFL+l_0)/betaT/l_0 + e*e*e*betaVL/betaU/l_0/g_0/g_0
    
    #print betaFL, g_1/g_0, h/l_0
    
    sigmaH = M/lamb/g_1/g_1
    sigmaR = (1.333*e*betaFL+l_0)*M/lamb/e/e/l_0
    sigmaTangenc = betaY*M/e/e - sigmaR*(K*K+1)/(K*K-1)
    sigma = {'sigmaR':sigmaR,
             'sigmaTangenc':sigmaTangenc,
             'sigmaH':sigmaH}
    return sigma
        
        
def mezenapeti(B,f,f_H):
    u'''
    Created on 17. 3. 2014
    11.5.4.2 - Fce vrátí hodnoty mezí jednotlivých napěti
    
    B - vnitřní průměr příruby
    f - dovolené namáhání
    f_H - dovolené namáhí krku ???
    @author: Stanislav Schütz
    '''
    if B <= 1000:
        k = 1
    elif B > 2000:
        k = 1.333
    else:
        k = 2.0/3*(1+B/2000)
        
    # Maximální napětí:
    sigmaH = 1.5*min(f,f_H)/k
    sigmaR = f/k
    # sigmaTangenc - hodnota je stejná jako sigmaR
    sigmaHplusR = 2*f/k
    # sigmaHplusTangenc - hodnota je stejná jako sigmaHplusR
    
    sigmamax = {'sigmaRmax':sigmaR,
             'sigmaHmax':sigmaH,
             'sigmaHplusR':sigmaHplusR}
    return sigmamax
        
        
def Vypocet(typ,geom,tesneni,materialVstup,zatizeni):
    # těsnění
    w = tesneni['w']
    y = tesneni['y']
    m = tesneni['m']
    
    # geometrie
    A = geom['A']
    B = geom['B']
    C = geom['C']
    G = geom['G']
    e = geom['e']
    g_0 = geom['g_0']
    g_1 = geom['g_1']
    h = geom['h']
    n = geom['n']
    d = geom['d']
    
    # zatížení
    P = zatizeni['P']
    T_prov = zatizeni['T_prov']
    T_mont = zatizeni['T_mont']
    
    # materiál
    mi = 0.3
    t = geom['t']
    matPriruba = materialVstup['matPriruba']
    
    # Vypočtení mat charakteristik pro dva stavy příruby
    provPriruba = Mat(matPriruba,t,T_prov)
    montPriruba = Mat(matPriruba,t,T_mont)
    
    # Přiřazení dovoleného namáhání k oběma stavům příruby
    f_prov = provPriruba['f']
    f_mont = montPriruba['f']
    
    if typ != 1: # dle typu příruby existuje i f_H pláště
        matPlast = materialVstup['matPlast'] 
        provPlast = Mat(matPlast,t,T_prov)
        montPlast = Mat(matPlast,t,T_mont)
        f_Hprov = provPlast['f']
        f_Hmont = montPlast['f']
    else:
        f_Hprov = f_prov
        f_Hmont = f_mont
        
    # podmínka pro materiál šroubů
    Re_BA = materialVstup['Re_BA']
    Re_B = materialVstup['Re_B']
    Rm = materialVstup['Rm']
    
    f_BA = min(Re_BA/3,Rm/4)
    f_B = min(Re_B/3,Rm/4)
    
    
    import math
    # kontrola šroubů
    zatizeni = Srouby(typ,w,n,d,g_1,B,C,G,P,f_BA,f_B,y,m)
    
    # Koeficienty pro výpočet
    deltaB = 2*C*math.sin(2*math.pi/n/2) # mnou spočtená vzdálenost mezi dvěmi šrouby
    C_F = max(math.sqrt(deltaB/(2*d + 6*e/(m+0.5))),1.0) # korekční součinitel rozteče šroubů
    K = A/B
    l_0 = math.sqrt(B*g_0)
    betaY = 1.0/(K-1)*(0.66845+5.7159*(K*K*math.log10(K))/(K*K-1.0))
    betaT = (K*K*(1 +8.55246*math.log10(K))-1)/(1.0472+1.9448*K*K)/(K-1)
    betaU = (K*K*(1 +8.55246*math.log10(K))-1)/(K*K-1)/1.36136/(K-1)
    
    
    def napeti(typ,B,g_0,g_1,e,h,K,l_0,betaT,betaU,betaY,mi,f,f_H,M,P,T):
        u'''Fce pro výpočet napětí v závislosti na momentu - řešení dvou stavů příruby'''
        
        # Rozhodnutí o způsobu výpočtu      
        if typ == 3:
            sigma = skrkem(g_0,g_1,e,h,K,l_0,betaT,betaU,betaY,mi,M) # volná metoda pro násuvnou přírubu s krkem
        elif g_0 <= 16 and P <= 2 and B/g_0 <= 300 and T <= 370:
            sigma = volna(e,betaY,M) # volná metoda pro typy 1 a 2
        else:
            sigma = integralne(g_0,g_1,e,h,K,l_0,betaT,betaU,betaY,mi,M) # integralni metoda
            
        sigmamax = mezenapeti(B,f,f_H) # získám maximální hodnoty napětí
        
        # Zapíšu hodnoty napětí a mezí napětí do výstupního slovníku
        napeti = sigmamax
        for clen in sigma:
            napeti[clen] = sigma[clen]
        
        return napeti

    
    if T_prov == T_mont:
        # Pro stejné teploty jsou stejné hodnoty f_prov a f_mont a tak stačí vybrat větší z momentů a vypočíst napětí jen jednou
        mom = max(zatizeni['M_op'],zatizeni['M_A'])
        M = mom*C_F/B
        vysl_prov = napeti(typ,B,g_0,g_1,e,h,K,l_0,betaT,betaU,betaY,mi,f_prov,f_Hprov,M,P,T_prov)
        vysl_mont = {}
    else:
        M = zatizeni['M_op']*C_F/B
        vysl_prov = napeti(typ,B,g_0,g_1,e,h,K,l_0,betaT,betaU,betaY,mi,f_prov,f_Hprov,M,P,T_prov)
        M = zatizeni['M_A']*C_F/B
        vysl_mont = napeti(typ,B,g_0,g_1,e,h,K,l_0,betaT,betaU,betaY,mi,f_mont,f_Hmont,M,P,T_prov)
    
    # Zápis všech výstupních napětí do slovníku
    vystup = {}
    for clen in zatizeni:
        vystup[clen] = zatizeni[clen]
    for clen in vysl_prov:
        vystup[clen+'_prov'] = vysl_prov[clen] # Index prov pro provozní hodnoty napětí
    for clen in vysl_mont:
        vystup[clen+'_mont'] = vysl_mont[clen] # A index mont pro napětí při montáži
    vystup['typ'] = typ
    vystup['n'] = n
    vystup['Re_prirMont'] = montPriruba['Re']
    vystup['Re_prirProv'] = provPriruba['Re'] # Do výpisu pouze Re
    if typ == 2 or typ == 3:
        vystup['Re_plastMont'] = montPlast['Re']
        vystup['Re_plastProv'] = provPlast['Re']
    return vystup