sher_duamel.py 7.66 KB
'''
Created on 28 авг. 2019 г.

@author: GadGet Life
'''

import math as m
import numpy as np
from matplotlib import mlab
import pylab
import csv
#import urllib.request
import requests
import json
import matplotlib.ticker as ticker
from fuzzylab import *
from datetime import datetime

# http://map.emercit.com/mgraph.php?mode=ultrasound&id=377&a=2020-01-10&b=2020-01-14&nocache=1578985253

class DuamelGidr:

    def __init__(self, time,amp,Tper):
        self.tbeg = time
        self.amp=amp
        self.T=Tper
        #print("U",time,amp)

    def calculation(self,tcur):

        t=tcur-self.tbeg
        return self.amp*(1-m.exp(-t/self.T))

def getLevel(q):
    Q=np.array([
        [0,    0],
        [0.345148257,    0.1],
        [1.097028578,    0.2],
        [2.158903614,    0.3],
        [3.491620674,    0.4],
        [5.071512601,    0.5],
        [6.882239402,    0.6],
        [8.911682045,    0.7],
        [11.15043894,    0.8],
        [13.59098985,    0.9],
        [16.22718582,    1],
        [19.05391666,    1.1],
        [22.06688277,    1.2],
        [25.26243236,    1.3],
        [28.63744147,    1.4],
        [32.18922318,    1.5],
        [35.91545733,    1.6],
        [39.81413517,    1.7],
        [43.88351498,    1.8],
        [48.12208599,    1.9],
        [52.52853871,    2],
        [57.10174025,    2.1],
        [61.84071361,    2.2],
        [66.74462006,    2.3],
        [71.81274423,    2.4],
        [77.04448112,    2.5],
        [82.43932505,    2.6],
        [87.99685986,    2.7],
        [93.71675045,    2.8],
        [99.59873533,    2.9],
        [105.64262,    3],
        [111.8482711,    3.1],
        [118.2156111,    3.2],
        [124.7446139,    3.3],
        [131.4353005,    3.4],
        [138.2877352,    3.5],
        [145.3020222,    3.6],
        [152.4783028,    3.7],
        [159.8167523,    3.8],
        [167.3175775,    3.9],
        [174.981015,    4]]
    )
    n=len(Q)
    for i in range(n):
        if q <= Q[i][0]:
            break
    return Q[i-1][1]+0.1*(q-Q[i-1][0])/(Q[i][0]-Q[i-1][0])    


class TestSherDuamel:

    def __init__(self, _k_up=2.8, _k_down=17.0, _k_multy=160.0, _k_addit=60.0, _zero_level=2.20):
        self.k_up = _k_up
        self.k_down = _k_down
        self.k_multy = _k_multy
        self.k_addit = _k_addit
        self.zero_level = _zero_level

        self.prec_old = 0.0
        self.prec_sum_old = 0.0

        self.UpperStore = []
        self.DownStore = []

    def CalculateLevel(self, _prec, _counter):
        prec_delta = _prec - self.prec_old

        print("prec_delta = {}".format(prec_delta))

        self.prec_old = _prec
        if np.abs(prec_delta) > 0.001:
            self.UpperStore.append(DuamelGidr(_counter, prec_delta, self.k_up))

        prec_sum = 0.0

        for UpStep in self.UpperStore:
            prec_sum = prec_sum + UpStep.calculation(_counter)

        print("prec_sum = {}".format(prec_sum))

        prec_sum_delta = prec_sum - self.prec_sum_old
        self.prec_sum_old = prec_sum

        if np.abs(prec_sum_delta) > 0.001:
            self.DownStore.append(DuamelGidr(_counter, prec_sum_delta, self.k_down))

        result_flow = 0.0

        for DownStep in self.DownStore:
            result_flow=result_flow+DownStep.calculation(_counter)

        print("result_flow = {}".format(result_flow))

        result_flow = result_flow * self.k_multy

        result_level = getLevel(result_flow + self.k_addit)-self.zero_level

        print("result_level = {}".format(result_level))

        return result_level




D=[]
M286=[]
normSUp=[]
normSDown=[]

L=[]
LUS=[]
count=0
with open('rs.csv', newline='') as csvfile:
    spamreader = csv.reader(csvfile, delimiter=';', quotechar='|')
    for row in spamreader:
        print(row[0])
        count=count+1
        D.append(datetime.strptime(row[0],'%d.%m.%y %H:%M'))
        M286.append(float(row[1]))
        L.append(float(row[4]))
        normSUp.append(float(row[2]))
        normSDown.append(float(row[3]))
        LUS.append(float(row[5]))

n_ya=0
o_ya=0

Tup=2.8;
Tdown=17.;

oldM286=0.
oldSUp=0.
oldSDown=0.

sumOoldDown=0.
sumOoldUp=0.
sumOold=0.0
Kkor=1.0
Kkor2=40.0
Kn=1
Tdob=9

Tcorr=6
Lpor=0.15
correction=True
k_filter=1.07
iosr=2

sdno=-0.7

DuamM286=[]
A=[]
B=[]
C=[]
D=[]
O=[]
U=[]
E=[]
F=[]

ND = []

dHold=0

shift=0
count=count-shift

test_duam = TestSherDuamel()

for i in range(count):

    # Получаем дельту в осадках на агк-286
    deltaM286=M286[i]-oldM286
    oldM286=M286[i]

    # Заполняем дюамелями с подъемными коэфами всё для осадков расчета по агк-286
    if np.abs(deltaM286) > 0.001:
        A.append(DuamelGidr(i,deltaM286,Tup))
        print(oldM286)
    sumO=0.0
    kf=len(A)

    for j in range(kf):
        sumO=sumO+A[j].calculation(i)
    sumO=Kkor*sumO
    deltaO=sumO-sumOold
    sumOold=sumO

    if np.abs(deltaO) > 0.001:
        B.append(DuamelGidr(i,deltaO,Tdown))

    O.append(sumO)

####################################################

    deltaM286=normSUp[i+shift]-oldSUp

    print("delta = {}".format(deltaM286))

    oldSUp=normSUp[i]

    if np.abs(deltaM286) > 0.001:
        C.append(DuamelGidr(i,deltaM286,Tup))
        print(oldM286)
    sumO=0.0
    kf=len(C)
    for j in range(kf):
        sumO=sumO+C[j].calculation(i)

    print("sum0 = {}".format(sumO))

    deltaO=sumO-sumOoldUp
    sumOoldUp=sumO

    print("delta0 = {}".format(deltaO))

    if np.abs(deltaO) > 0.001:
        D.append(DuamelGidr(i,deltaO,Tdown))

    sumB=0.0
    kf=len(D)
    for j in range(kf):
        sumB=sumB+D[j].calculation(i)

    # test

    print("sumB = {}".format(sumB))

    sumB=sumB*160.0

    sumB=getLevel(sumB+60)-2.20

    print("LEVEL= = {}".format(sumB))

    U.append(sumB)

    test_result = test_duam.CalculateLevel(normSUp[i+shift], i)
    ND.append(test_result)

    if sumB != test_result:
        print("gopa!")
        print("gopa!")

    print("old:{}, new:{}, RESULT={}".format(sumB, test_result, sumB==test_result))

####################################################

    deltaM286=normSDown[i]-oldSDown
    oldSDown=normSDown[i]
    if np.abs(deltaM286) > 0.001:
        E.append(DuamelGidr(i,deltaM286,Tup))
    sumO=0.0
    kf=len(E)
    for j in range(kf):
        sumO=sumO+E[j].calculation(i)
    sumO=sumO


    deltaO=sumO-sumOoldDown
    sumOoldDown=sumO

    if np.abs(deltaO) > 0.001:
        F.append(DuamelGidr(i,deltaO,Tdown))

xlist = np.array([range(0,10*count,10)]).T


fig, ax = pylab.subplots()

ax.plot (xlist, U,label="RadarUp", linewidth = 2)
ax.plot (xlist, ND,label="New_Duamel", linewidth = 2)
ax.plot (xlist, L, label="Real Level", linewidth = 3)


ax.minorticks_on()

ax.xaxis.set_major_locator(ticker.MultipleLocator(300))
ax.xaxis.set_minor_locator(ticker.MultipleLocator(60))

ax.grid(which='major',
        linewidth = 1)
ax.grid(which='minor',linestyle = ':')

ax.legend()

#ax.hlines(n_ya, 0, 10*count, color = 'y', linestyle = '--')
#ax.hlines(o_ya, 0, 10*count, color = 'r', linestyle = '--')
#ax.hlines(-4, -5, 5)
ax.set_xlabel('Time (minute)')
ax.set_ylabel('Level (m)')

#st="Tup:{}   Tdown:{}    date:{}-{}    Correction:{}    Tdob:{}    Kmain:{}   Lpor:{}"

#ax.set_title(st.format(str(Tup),str(Tdown),dbeg,dend,correction,Tdob,Kkor,Lpor))



#fig.set_figwidth(18)
#fig.set_figheight(12)

pylab.show()
#sf="out_{}_{}_{}_{}.pdf"
#fig.savefig(sf.format(dbeg,dend,str(Tup),str(Tdown)))