<jupyter_start><jupyter_text><jupyter_code>from google.colab import files
data_to_load = files.upload()<jupyter_output><empty_output><jupyter_text>New Section<jupyter_code>import xlrd 
import numpy as np  
# Give the location of the file 
loc = ("tag1.xlsx") 
# To open Workbook 
wb = xlrd.open_workbook(loc) 
sheet = wb.sheet_by_index(0) 
# For row 0 and column 0 
ai=sheet.cell_value(1, 0)
Numrow=sheet.nrows
fi=1
for  i in range(Numrow):
  ain=sheet.cell_value(i, 0)
  if abs(ain-ai)>0.00001:
     fi+=1
     ai=ain   
fj=int(np.round(Numrow/fi))
matrovtp=np.zeros((fi,fj))
matvis=np.zeros((fi,fj))
chv=0  
brt=np.zeros((1,fj))
brp=np.zeros((1,fi))
for i in range(fi):
  brp[0,i]=sheet.cell_value(chv, 0)
  for j in range(fj):
    matrovtp[i,j]=sheet.cell_value(chv, 2)
    matvis[i,j]=sheet.cell_value(chv,3)    
    if i==0:
      brt[0,j]=sheet.cell_value(chv, 1)
    chv+=1  
brp
p0ta=sheet.cell_value(0, 0)
t0ta=sheet.cell_value(0, 1)
dpta=sheet.cell_value(fi, 0)-sheet.cell_value(fi-1, 0)
dtta=sheet.cell_value(1, 1)-sheet.cell_value(0, 1)
Taz=40
paz=10000000
Npta=int((paz-p0ta)/dpta)
Ntta=int((Taz-t0ta)/dtta)
rowp=matrovtp[Npta,Ntta]+(paz-brp[0,Npta])*(matrovtp[Npta+1,Ntta]-matrovtp[Npta,Ntta])/dpta+(Taz-brt[0,Ntta])*(matrovtp[Npta,Ntta+1]-matrovtp[Npta,Ntta])/dtta
from pylab import *        
#plt.plot(matxl[:,0],'g')
plt.plot( brp.T/(matrovtp[:,49:50]*437*(brt[0,49]+273.3)),'r')
from pylab import *        
#plt.plot(matxl[:,0],'g')
###plt.plot( brp.T/(matrovtp[:,49:50]*437*(brt[0,49]+273.3)),'g')
#plt.plot(matvis[1,10],'g')
import xlrd 
import numpy as np  
# Give the location of the file 
#loc = ("pcc1000.xlsx") 
#loc = ("pt1.xlsx") #q=100
#loc = ("PT50,1000.xlsx") 
loc = ("PTl150.1000.xlsx")
#loc = ("PT80.1000.xlsx")
#loc =("PTL100.110.xlsx")
#loc = ("PTW100.xlsx") #q=100
#loc = ("ptl000.xlsx") #q=100
#loc = ("PTS1.xlsx") #q=100
# To open Workbook 
#loc = ("uqc1000.xlsx") 
wb = xlrd.open_workbook(loc) 
sheet = wb.sheet_by_index(0) 
# For row 0 and column 0 
#ai=sheet.cell_value(0, 0)
Numrow=sheet.nrows
LPT=np.zeros((Numrow,3))
NX=50
Tempo=np.zeros((NX+1,1))
L=50000
nol=L/(Numrow-1)
for  i in range(Numrow):
  LPT[i,0]=sheet.cell_value(i, 0)
  LPT[i,1]=sheet.cell_value(i, 1)
  LPT[i,2]=sheet.cell_value(i, 2)  
import xlrd 
import numpy as np  
# Give the location of the file 
#loc = ("uqc1000.xlsx") 
loc = ("QUL150-1000.xlsx") 
#loc = ("QUS (1).xlsx") 
#loc = ("QUl1000.110.xlsx") 
#loc = ("QU80.1000.xlsx") 
# To open Workbook 
wb = xlrd.open_workbook(loc) 
sheet = wb.sheet_by_index(0) 
# For row 0 and column 0 
ai=sheet.cell_value(0, 0)
Numrow=sheet.nrows
LUQ=np.zeros((Numrow,3))
NX=50
Tempo=np.zeros((NX+1,1))
L=50000
nol=L/(Numrow-1)
for  i in range(Numrow-1):
  LUQ[i,0]=sheet.cell_value(i, 0)
  LUQ[i,1]=sheet.cell_value(i, 1)
  LUQ[i,2]=sheet.cell_value(i, 2)

#tart100Nas0
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 10 17:33:55 2020
@author: Saeed
"""
import numpy as np
#import random as ran
import math as ma
NXT=160
NX=NXT
NT=300
#NT=1000
disx=LPT[180,0]-LPT[20,0]
dx=disx/NX
prdl=1
hj=np.zeros((prdl,1))
pp=np.zeros((prdl,2))
ptot=np.zeros((NX+1,1))
rovhol=np.zeros((NX+1,1))
Uhold=np.zeros((NX,1))
matxl=np.zeros((prdl,1))
beta=0.1
attopa=101.32533538686013
CRPRE=66.1799*101325.33538686013
CRTEMP=214.846
az=0.257/CRPRE
bz=0.533*(CRTEMP/CRPRE)
p0=20802484
p0=7802484
p0s=LPT[20,1]*10**5
R=467
R=437
#R=520
#R=460 
upon=0
gh=0
absvis=1.2e-04
absvis=1e-05
rog=0.0008
D=0.75
A=(D/2)**2*ma.pi
uy=np.zeros((2,2721))
yui=0
g=9.8
vi=np.zeros((1,NX))
viv=np.zeros((1,NX))
#dx=disx/NX
def solfr(D,rog,Re):
  import math as ma
  te=0
  f2=0.19593676
  f1=f2
  for i in range(1):
     te=1/f1**0.5+2*ma.log(rog/(3.7*D)+2.51/(Re*f1**0.5),10)
     f1=f1+0.01*te
     f1=0.25/(ma.log(rog/(3.7*D)+5.74/((Re)**0.9),10))**2
     if f1<0:
       f1=f2/10
     if te<0.00001:
        break  
  return f1        
#mdatNl=1000 
tempe=LPT[180,2]
temps=LPT[20,2] 
shibmat=np.zeros((NXT,1))
ypipe=np.zeros((NXT+1,1))
yp=np.zeros((int(NXT*dx)+1,1))
#shibmat= 0.1*(np.random.random((NX,1))-0.5)
for i in range(1,NX+1):
  ypipe[i,0]=ypipe[i-1,0]+shibmat[i-1,0]*dx
  for j in range(int(dx)):
     yp[j+(i-1)*int(dx),0]=j*((ypipe[i,0]-ypipe[i-1,0])/dx)+ypipe[i-1,0]
 #  xl=disx/2
xl=disx
 #NX=int(np.round((disx-xl)*NXT/(disx)))+1
NX=NX
##NX=48
NX2=NX
mdat=1000     
pout=5015727
pout=LPT[180,1]*10**5
Temp1=tempe*(NX-NX)/NXT+(NX/NXT)*temps  #plus 273.73
Temp2=tempe
#disx1=disx-xl
disx1=xl
#dx=disx1/NX
landa=1000
P=np.zeros((NX+1,NT+1))
U=np.zeros((NX+1,NT+1))
pmat=np.zeros((NX-1,NX-1))
rov=np.zeros((NX+1,NT+1))
rovm=np.zeros((NX,NT+1))
Temp=np.zeros((NX+1,1))
z=np.zeros((NX+1,1))
#vc=np.zeros((NT,1))
p0mat=np.zeros((1,NT))
b=np.zeros((NX-1,1))
bsv=np.zeros((NX,1))
Pc=np.zeros((NX,1))
d=np.zeros((NX,1))
mdat1=np.zeros((NX,NT))
lap=1
#p0=5010727
#p0=90102484
yui=0
inrov=p0/(R*312.485)
endrov=pout/(R*312.485)    
p01=p0
landa=1000
Tempp=np.zeros((NX+1,1))
for  i  in range (NX+1):
      Temp[i,0:1]=i*((Temp2-Temp1)/NX)+Temp1
rovm0=rovm[0,0:1]
ppi=0 
dif2=np.zeros((NX,1))
j=0
Cv=1.85
Cp=2.34
kc=Cp/Cv
amatr=np.zeros((NX,NX))
bmat=np.zeros((1,NX))
alfac=0
npi=np.pi
ohtc=8
temamb=50
Zc=np.zeros((NX+1,1))
Zc12=np.zeros((NX+1,1))
Umin=np.zeros((NX,NT))
e0=np.zeros((NT,1))
asd=0
asp=0
alfac=1
alfap=1
ind=1
mattfr=np.zeros((1,NX))
pf=np.zeros((1,NX))
pout=LPT[180,1]*10**5
P[NX,0]=pout
Temp[NX,0:1]=LPT[180,2]
Tempp[NX,0:1]=LPT[180,2]
vi=np.zeros((1,NX+1)) 
tinf=8
ham=8
mcs=50    #W/(m.K) materials conductivity steel
mcc=1    #W/(m.K)  materials conductivity concrete
Lst=0.009  #lengh steel
Lco=0.02   #lengh concrete
chtc=8 #W/(m**2.K)  convection heat transfer coeficient 
ohtc= 1/((1/8)+(Lco/1)+(Lco/1)+(Lst/50)) #W/(m**2.K)   overal heat transfer coeficient 
Rin=1/(ham*pi*D*dx)
Rst=np.log((D+2*Lst)/D)/(2*pi*mcs*dx)
Rcon1=np.log((D+2*Lst+2*Lco)/(D+2*Lst))/(2*pi*mcs*dx)
Rcon2=np.log((D+2*Lst+2*Lco+2*Lco)/(D+2*Lst+2*Lco))/(2*pi*mcs*dx)
Rair=1/(pi*(D+2*Lst+2*2*Lco)*dx*10)
Rtot=Rin+Rst+Rcon1+Rcon2+Rair
ps=LPT[20,1]*10**5
asd1=0
for j in range(0,NT):  
             P[NX,j]=pout
             for iv in range(NX,-1,-1):  
                if j>1:          
                    Temp[1:NX,0]=0.8*Tempp[1:NX,0]+0.2*Temp[1:NX,0]                                                      
                Taz=Temp[iv,0]
                paz=P[iv,j]
                Npta=51
                Ntta=51
                for ip in range(fi-1):
                  if brp[0,ip+1]>paz>=brp[0,ip]:
                    dpta=brp[0,ip+1]-brp[0,ip]
                    Npta=ip       
                    break
                for it in range(fj-1):     
                  if brt[0,it+1]>Taz>=brt[0,it]:
                    dtta=brt[0,it+1]-brt[0,it]
                    Ntta=it          
                    break
                if (Npta>fi or Ntta>fj): 
                      rov[iv,j]=(P[iv,j])/(R*(Temp[iv,0]+273.3)) 
                      absvis=0.00005
                else:
                     odif=(paz-brp[0,Npta])*(matrovtp[Npta+1,Ntta]-matrovtp[Npta,Ntta])/dpta+(Taz-brt[0,Ntta])*(matrovtp[Npta,Ntta+1]-matrovtp[Npta,Ntta])/dtta
                     ttdif=0.5*(Taz-brt[0,Ntta])**2*(matrovtp[Npta,Ntta+1]-2*matrovtp[Npta,Ntta]+matrovtp[Npta,Ntta-1])/dtta**2
                     pdif=0.5*(paz-brp[0,Npta])**2*(matrovtp[Npta+1,Ntta]-2*matrovtp[Npta,Ntta]+matrovtp[Npta-1,Ntta])/dpta**2
                     rov[iv,j]=matrovtp[Npta,Ntta]+pdif+ttdif+odif
                     onedif=(paz-brp[0,Npta])*(matvis[Npta+1,Ntta]-matvis[Npta,Ntta])/dpta+(Taz-brt[0,Ntta])*(matvis[Npta,Ntta+1]-matvis[Npta,Ntta])/dtta
                     absvis=matvis[Npta,Ntta]+onedif         
                Re=(mdat*D)/(A*absvis) 
                landa=0.25/(np.log10((rog/(3.7*D))+(5.74/(Re**0.9))))**2
                if iv>0 and   iv<NX-1 and j>5 :
                   dif2=-2*U[iv,j-1]+U[iv-1,j-1]+U[iv+1,j-1]
                   ind=(U[iv-1,j-1]+(1/12)*dif2)/ U[iv-1,j-1]
                fr=ind*landa*(mdat)*(dx)/(2*D)         
                U[iv,j]=mdat/(A*rov[iv,j])
                if j==0:
                   if iv>0:
                      P[iv-1,j]=P[iv,j]+fr*U[iv,j]/A        
                else:
                   if iv>0:
                      Umin[iv-1,j]=(U[iv,j-1]+U[iv-1,j-1])/2
                   Zc[:,0]=P[:,j-1]/(rov[:,j-1]*437*(Temp[:,0]+273.3))
                   if iv>0:
                      P[iv-1,j]=P[iv,j]+fr*Umin[iv-1,j]/A+1*mdat*(U[iv,j-1]-U[iv-1,j-1])/A+alfap*(mdat*(U[iv,j-1]-U[iv-1,j-1])/A)
                      vi[0,iv-1]=alfac*landa*((dx)*Umin[iv-1,j]**2)/(Cp*1000*D)
                     #  vi[0,iv-1]=0.5*Umin[iv-1,j]*(U[iv,j]-U[iv-1,j])/(Cp*1000)
                      tmin=(Tempp[iv-1,0]+Tempp[iv,0])/2
                      rovmin=(rov[iv-1,0]+rov[iv,0])/2
                      tl=dx/(Umin[iv-1,j]) 
                      #Tempp[iv-1,0]=(Tempp[iv,0]+273.15)*(P[iv-1,j]/P[iv,j])**((kc-1)/kc)-273.15-vi[0,iv-1]-0.01*D*pi*dx*(tmin-30)/(1000*Cp*Rtot*mdat)
                      Tempp[iv-1,0]=(Tempp[iv,0]+273.15)*(P[iv-1,j]/P[iv,j])**((kc-1)/kc)-273.15-vi[0,iv-1]
                      # Tempp[iv-1,0]=LPT[20+iv-1,2]
                      
                   else:  
                      asd=(0.001*asd+0.005*(Tempp[0,0]-temps))          
                      alfac=asd+alfac 
                      asd1=0.1*asd1+(100/ps)*(P[0,j-1]-ps)          
                      alfap=-asd1+alfap 
                      alfap=0 
                          
from pylab import *        
#plt.plot(matxl[:,0],'g')
plt.plot(P[:,j],'g')
#plt.plot(ypipe[:,0],'g')
alfac


#plt.plot(LPT[20:181,2],'g')
#plt.plot(LPT[20:181,2]-Tempp[:,0],'g')
#plt.plot(P[:,j],'r')
#plt.plot(Tempp[:,0],'r')
plt.plot(-P[:,j]+LPT[20:181,1]*10**5,'g')
#plt.plot(fgh,'g')
#fgh=P[:,j]-LPT[20:181,1]*10**5
#plt.plot(,'r')
#plt.plot(thh[:,0],'r')
#plt.plot(P[:,j],'r')
#plt.plot(P[:,j],'r')
#LPT[180,2]
#(61.8302917480469+25.6020698547363)/2
#P[0,j-1]
#tafantz
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 10 17:33:55 2020
@author: Saeed
"""
import numpy as np
#import random as ran
import math as ma
NXT=160
NT=100
prdl=100
fi=50
fj=50
ptot=np.zeros((NXT+1,1))
rovhol=np.zeros((NXT+1,1))
Uhold=np.zeros((NXT+1,1))
Thold=np.zeros((NXT+1,1))
matxl=np.zeros((prdl,1))
disx=LPT[180,0]-LPT[20,0]
dx=disx/NX
p0=23021755.981445347
R=460 
rog=0.0008
D=0.75
A=(D/2)**2*ma.pi
mdatNl=LUQ[180,1]
#mdatNl=776.2
mdatbl=LUQ[20,1]
g=9.8
dx=disx/NXT
NX=NXT
vi=np.zeros((1,NXT+1))
p0ta=1000000
t0ta=0
dpta=brp[0,1]-brp[0,0]
dtta=brt[0,1]-brt[0,0]
Cv=1.85
Cp=2.34
kc=Cp/Cv
p0s=LPT[20,1]*10**5
p0=LPT[20,1]*10**5
tempe=0
def solfr(D,rog,Re):
  import math as ma
  te=0
  f2=0.19593676
  f1=f2
  for i in range(1):
     te=1/f1**0.5+2*ma.log(rog/(3.7*D)+2.51/(Re*f1**0.5),10)
     f1=f1+0.01*te
     f1=0.25/(ma.log(rog/(3.7*D)+5.74/((Re)**0.9),10))**2
     if f1<0:
       f1=f2/10
     if te<0.00001:
        break  
  return f1        
#landabl=solfr(D,rog,Rebl)
#landaNl=solfr(D,rog,ReNl)
tempe=LPT[180,2]
temps=LPT[20,2]
shibmat=np.zeros((NXT,1))
ypipe=np.zeros((NXT+1,1))
yp=np.zeros((int(NXT*dx)+1,1))
alfac=0.5385624101071541
pp=zeros((prdl,2))
TT=zeros((prdl,2))
bnj=np.zeros((1,prdl))
pout=5231451.035
ind=1
yui=0
#shibmat= 0.1*(np.random.random((NXT,1))-0.5)
alfac2=0.0003
alfac1=0.317
tinf=8
ham=8
mcs=50    #W/(m.K) materials conductivity steel
mcc=1    #W/(m.K)  materials conductivity concrete
Lst=0.009  #lengh steel
Lco=0.02   #lengh concrete
chtc=8 #W/(m**2.K)  convection heat transfer coeficient 
ohtc= 1/((1/8)+(Lco/1)+(Lco/1)+(Lst/50)) #W/(m**2.K)   overal heat transfer coeficient 
Rin=1/(ham*pi*D*dx)
Rst=np.log((D+2*Lst)/D)/(2*pi*mcs*dx)
Rcon1=np.log((D+2*Lst+2*Lco)/(D+2*Lst))/(2*pi*mcs*dx)
Rcon2=np.log((D+2*Lst+2*Lco+2*Lco)/(D+2*Lst+2*Lco))/(2*pi*mcs*dx)
Rair=1/(pi*(D+2*Lst+2*2*Lco)*dx*10)
Rtot=Rin+Rst+Rcon1+Rcon2+Rair
hyu=0

for i in range(1,NXT+1):
  ypipe[i,0]=ypipe[i-1,0]+shibmat[i-1,0]*dx
  for j in range(int(dx)):
     yp[j+(i-1)*int(dx),0]=j*((ypipe[i,0]-ypipe[i-1,0])/dx)+ypipe[i-1,0]   
for l in range(prdl):
    asd=0
    asd2=0
    if l==0:
        xl=disx/2
       #  xl=disx
        grl=0
    else:
       grl=(disx/(2*p0))*(p0-pp[l-1,1])+0.5*grl
       xl=xl+grl
    matxl[l,0]=xl   
    for zl in range(2):        
        if zl==0: 
            NX=int(np.round((disx-xl)*NXT/(disx)))+1
            #NX=50
            ##NX=48
            NX2=NX
            mdat=mdatNl             
            pout=LPT[180,1]*10**5
            tout=tempe
            disx1=disx-xl
             #disx1=xl
            dx=disx1/NX
          #  landa=landaNl
            Temp1=tempe*(NXT-NX)/NXT+(NX/NXT)*temps  #plus 273.73
            Temp2=tempe
        else:
            disx1=xl
            Temp2=Temp[0,0]
            NX=NXT-NX2
            mdat=mdatbl
            # pout= pp[l,0]+(Uhold[NX-1,0]-Uhold[NX,0])*mdatNl/A+ Uhold[NX-1,0]*mdatbl/A
            pout= pp[l,0] 
            tout=TT[l,0]    
       #     landa=landabl       
            dx=xl/NX  
            Temp2=Temp[0,0]
            Temp1=temps   
        P=np.zeros((NX+1,NT+1))
        U=np.zeros((NX+1,NT+1))

        rov=np.zeros((NX+1,NT+1))
        rovm=np.zeros((NX,NT+1))
        Umin=np.zeros((NX,NT+1))
        Temp=np.zeros((NX+1,1))
        Tempp=np.zeros((NX+1,1))
        z=np.zeros((NX+1,1))
        #vc=np.zeros((NT,1))
        p0mat=np.zeros((1,NT))
        b=np.zeros((NX-1,1))
        bsv=np.zeros((NX,1))
        Pc=np.zeros((NX,1))
        d=np.zeros((NX,1))
        mdat1=np.zeros((NX,NT))
        Temp[NX,0]=Temp1
        Tempp[NX,0]=Temp1
   #     if l<1:        
   #         for  i  in range (NX+1):
   #             Temp[i,0:1]=i*((Temp2-Temp1)/NX)+Temp1
   #             Tempp[i,0:1]=i*((Temp2-Temp1)/NX)+Temp1

  #      else:
  #          if zl==0:
  #              Tempp[:,0]=Thold[NXT-NX:NXT+1,0]
  #              Temp[:,0]=Thold[NXT-NX:NXT+1,0]
                
  #          else:  
  #              Tempp[:,0]=Thold[0:NX+1,0]
  #              Temp[:,0]=Thold[NXT-NX:NXT+1,0]            
      #      for  rav  in range (NX+1):
      #          Temp[rav,0:1]=rav*((Temp2-Temp1)/NX)+Temp1  
        for j in range(0,NT):   
             Tempp[NX,0]=tout  
             if j<1:        
               for  i  in range (NX+1):
                  Temp[i,0:1]=i*((Temp2-Temp1)/NX)+Temp1
                  Tempp[i,0:1]=i*((Temp2-Temp1)/NX)+Temp1 
             P[NX,j]=pout
             
             for iv in range(NX,-1,-1):  
                if j>2:          
                    Temp[1:NX,0]=0.1*Tempp[1:NX,0]+0.9*Temp[1:NX,0]   
                    Temp=Tempp                                                                                                     
                Taz=Temp[iv,0]
                paz=P[iv,j]
                Npta=51
                Ntta=51
                for ip in range(fi-1):
                  if brp[0,ip+1]>paz>=brp[0,ip]:
                    dpta=brp[0,ip+1]-brp[0,ip]
                    Npta=ip       
                    break
                for it in range(fj-1):     
                  if brt[0,it+1]>Taz>=brt[0,it]:
                    dtta=brt[0,it+1]-brt[0,it]
                    Ntta=it          
                    break
                if (Npta>fi or Ntta>fj): 
                      rov[iv,j]=(P[iv,j])/(R*(Temp[iv,0]+273.3)) 
                      absvis=0.00005
                else:
                     odif=(paz-brp[0,Npta])*(matrovtp[Npta+1,Ntta]-matrovtp[Npta,Ntta])/dpta+(Taz-brt[0,Ntta])*(matrovtp[Npta,Ntta+1]-matrovtp[Npta,Ntta])/dtta
                     ttdif=0.5*(Taz-brt[0,Ntta])**2*(matrovtp[Npta,Ntta+1]-2*matrovtp[Npta,Ntta]+matrovtp[Npta,Ntta-1])/dtta**2
                     pdif=0.5*(paz-brp[0,Npta])**2*(matrovtp[Npta+1,Ntta]-2*matrovtp[Npta,Ntta]+matrovtp[Npta-1,Ntta])/dpta**2
                     rov[iv,j]=matrovtp[Npta,Ntta]+pdif+ttdif+odif
                     onedif=(paz-brp[0,Npta])*(matvis[Npta+1,Ntta]-matvis[Npta,Ntta])/dpta+(Taz-brt[0,Ntta])*(matvis[Npta,Ntta+1]-matvis[Npta,Ntta])/dtta
                     absvis=matvis[Npta,Ntta]+onedif
             
                Re=(mdat*D)/(A*absvis) 
                landa=solfr(D,rog,Re)
               # landa=0.25/(np.log10((rog/(3.7*D))+(5.74/(Re**0.9))))**2
                if iv>0 and   iv<NX-1 and j>2 :
                   dif2=-2*U[iv,j-1]+U[iv-1,j-1]+U[iv+1,j-1]
                   ind=(U[iv-1,j-1]+(1/12)*dif2)/ U[iv-1,j-1]
                fr=ind*landa*(mdat)*(dx)/(2*D)
                
                U[iv,j]=mdat/(A*rov[iv,j])
                if j==0 and l==0:
                   if iv>0:
                      P[iv-1,j]=P[iv,j]+fr*U[iv,j]/A 
                    
                else:
                   
                  # Zc[:,0]=P[:,j-1]/(rov[:,j-1]*437*(Temp[:,0]+273.3))
                   if iv>0:
                      Umin[iv-1,j]=(U[iv,j-1]+U[iv-1,j-1])/2
                      P[iv-1,j]=P[iv,j]+fr*Umin[iv-1,j]/A +mdat*(U[iv,j-1]-U[iv-1,j-1])/A-3.334118751019441*(mdat*(U[iv,j-1]-U[iv-1,j-1])/A)   
                      if zl==0:
                        #  alfac=alfac1
                          alfac=0.28639507516523577
                      else:
                       #  alfac=0.317 
                         alfac=0.28639507516523577
                      vi[0,iv-1]=alfac*landa*((dx)*Umin[iv-1,j]**2)/(Cp*1000*D)  
                      vi[0,iv-1]
                      #vi[0,iv-1]=alfac*(fr*Umin[iv-1,j])/((rov[iv-1,j]+rov[iv,j])*A*Cp*500)
                    
                  #    texp=np.exp((pi*D*dx*Umin[iv-1,j])/(1000*Cp*mdat))
                  #    tam=tinf-(tinf-Tempp[iv,0])*np.exp((pi*D*dx*ohtc)/(1000*Cp*mdat))-Tempp[iv,0]
                      tmin=(Temp[iv-1,0]+Temp[iv,0])/2
                      rovmin=(rov[iv-1,0]+rov[iv,0])/2
                      tl=dx/(Umin[iv-1,j])
                      if   zl==1: 
                            Tempp[iv,0]=Tempp[iv,0]+alfac2*iv
                      else: 
                            Tempp[iv,0]=Tempp[iv,0]+alfac2*(NX-iv) 
                      Tempp[iv-1,0]=(Tempp[iv,0]+273.15)*(P[iv-1,j]/P[iv,j])**((kc-1)/kc)-273.15-vi[0,iv-1]#(tmin-30)/(1000*Cp*Rtot*mdat)+
                     # if zl==0 and j>0:
                      #   Tempp[:,0]=Thold[NXT-NX:NXT+1,0]
                     # else:   
                      #    Tempp[:,0] = Thold[0:NX+1,0]
                   #   Tempp[iv-1,0]=(Tempp[iv,0]+273.15)*(P[iv-1,j]/P[iv,j])**((kc-1)/kc)-273.15-vi[0,iv-1]+tam
                    #  Tempp[iv-1,0]=LPT[iv+19,2]
                     # Tempp[iv-1,0]=(Tempp[iv,0]+273.15)
                   else:  
                     if zl==1 and j==NT-1 and l>6:
                        asd=(0.001*asd+0.001*(TT[l-1,zl]-temps))          
                        alfac1=alfac1+asd
                        asd2=(0.00001*asd2+0.00001*(TT[l-1,zl]-temps)) 
                        asd2=(0.00001*asd2+0.00001*(Tempp[0,0] -temps))
                              
                        alfac2=alfac2-asd2
                        hyu+=1
        if zl==0:
            ptot[NXT-NX:NXT+1,0]=P[:,j] 
            rovhol[NXT-NX:NXT+1,0]=rov[:,j]
            Uhold[NXT-NX:NXT+1,0]=U[:,j]
            Thold[NXT-NX:NXT+1,0]=Tempp[:,0]
        else:
            ptot[0:NX+1,0]=P[:,j]
            rovhol[0:NX+1,0]=rov[:,j]
            Uhold[0:NX+1,0]=U[:,j]   
            Thold[0:NX+1,0]=Tempp[:,0]   
     #   Thold[:,0]=LPT[20:181,2]         
        pp[l,zl]=P[0,j]
        TT[l,zl]=Tempp[0,0]     
from pylab import *        
plt.plot(matxl[:,0],'g')
#plt.plot(P[:,j],'g')
#plt.plot(ypipe[:,0],'g')
plt.ylabel('length,(km)')
plt.xlabel('number-repeat')
legend(['model'])
import matplotlib.pyplot as plt
twop=np.zeros((NX))
twop=range(0,NX)
Xj = np.arange(0.125, 40.125,0.25 )

Xj = np.arange(0, 40.250,0.25 )

twop=Xj
LPT[20:181,1]*10**5
####plt.plot(twop,Tempp[:,0],'b',linestyle='dashed')
#plt.plot(LUQ[21:181,0]/1000,LUQ[21:181,2],'b')
#plt.plot(LPT[20:181,0]/1000 ,LPT[20:181,1]*10**5,'b')
#plt.plot(LPT[20:181,0]/1000 ,LPT[20:181,2],'b')
#plt.plot(, Uhold[:,0],'r',linestyle='dashed')
#plt.plot(twop,LPT[20:181,2],'b')
#
##plt.plot(LPT[20:181,0]/1000,Thold[:,0],'r',linestyle='dashed')
#plt.plot(LPT[20:181,0]/1000,Uhold[:,0],'r',linestyle='dashed')
#plt.plot(twop,rovhol[:,0],'g')
#plt.plot(twop,LPT[20:181,1],'b')
#plt.plot(twop,[:,0],'r',linestyle='dashed')
#plt.plot(LPT[20:181,0]/1000,Tempp[:,0],'r',linestyle='dashed')
#plt.plot(LPT[20:181,0]/1000,U[:,j-1],'r',linestyle='dashed')
#plt.plot(Xj ,rovhol[:,0],'b')
#plt.plot(twop,Temp[:,0],'red')
#plt.plot(twop,ptot[:,0]-LPT[20:181,1]*10**5,'r',linestyle='dashed')
#plt.plot(twop,LPT[20:181,2]-Temp[:,0],'r')
#plt.plot(LPT[20:181,0]/1000,LPT[20:181,1]*10**5-P[:,j],'g')
#plt.plot(LPT[20:181,0]/1000,tre1,'b')
#plt.plot(LPT[20:181,0]/1000,LPT[20:181,2]-Temp[:,0],'g')
#plt.plot(twop,tre1,'r')
#tre1=LPT[20:181,1]*10**5-P[:,j]
#plt.plot(twop,Umin[:,j],'r',linestyle='dashed')
#plt.ylabel('GAS mass flow,(kg/s)')
#plt.ylabel('errore,(Pa)')
#plt.ylabel('errore,()')
#plt.ylabel('pressure,(Pa)')
#plt.ylabel('velocity,(m/s)')
plt.ylabel('temperature,(°C)')
plt.xlabel('length,(km)')
#legend(['Correction.fr','primitive.fr' ])
#legend(['OLGA-model'])
#legend(['OLGA','model'])
#legend(['E.Pcorrect','E.Pprime' ])
#plt.savefig("ftyh.jpg", dpi=1000)
#ftyh
plt.plot(matxl[:,0]+LPT[20,0],'g')
LPT[20,0]/1000<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2277596866298551, 'fraction_numerical': 0.09589881542309689, 'mean_word_length': 2.3047411104179663, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '428760', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 11291, 'n_tokens_neox': 10271, 'n_words': 1001, 'path': 'payanname.ipynb', 'repo_name': 'saeed13579/docs', 'size': '21189'}