text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: HeatExchangerDetailed.mso 197 2007-03-08 14:31:57Z bicca $
#*------------------------------------------------------------------
type ShellandTubesBasic
	ShellandTubesBasic()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Hot Side in the Exchanger",
				:Hidden=>true,
				:Valid=>["shell","tubes"],
				:Default=>"shell"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"TEMA Designation",
				:Valid=>["Eshell","Fshell"],
				:Default=>"Eshell"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Tube Layout Characteristic Angle",
				:Valid=>["Triangle","Rotated Square","Square"],
				:Default=>"Triangle"
			)),
			Tube_Side_Main(Dict{Symbol,Any}(
				:Brief=>"Tube Side Section" ,
				:Symbol=>"^{tube}"
			)),
			Shell_Side_Main(Dict{Symbol,Any}(
				:Brief=>"Shell Side Section" ,
				:Symbol=>"^{shell}"
			)),
			Baffles_Main(Dict{Symbol,Any}(
				:Brief=>"Baffle Section",
				:Symbol=>" "
			)),
			Clearances_Main(Dict{Symbol,Any}(
				:Brief=>"Diametral Clearances",
				:Symbol=>" "
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet Tube Stream",
				:PosX=>0.13,
				:PosY=>0,
				:Symbol=>"_{in }^{tube}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet Tube Stream",
				:PosX=>0.13,
				:PosY=>1,
				:Symbol=>"_{out }^{tube}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet Shell Stream",
				:PosX=>0.31,
				:PosY=>1,
				:Symbol=>"_{in }^{shell}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet Shell Stream",
				:PosX=>0.81,
				:PosY=>0,
				:Symbol=>"_{out }^{shell}"
			)),
			Details_Main(Dict{Symbol,Any}(
				:Brief=>"Details in Heat Exchanger",
				:Symbol=>" "
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Number of Tube rows Crossed in one Crossflow Section",
				:Hidden=>true,
				:Lower=>1
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Number of Effective Crossflow rows in Each Window",
				:Hidden=>true,
				:Lower=>1
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Variable for calculating Ji heat transfer correction Factor",
				:Hidden=>true,
				:Lower=>1e-3
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Variable for calculating shell side pressure drop friction Factor",
				:Hidden=>true,
				:Lower=>1e-3
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"ByPass Correction Factor for Pressure Drop",
				:Hidden=>true,
				:Lower=>1e-3
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Correction Factor for Pressure Drop",
				:Hidden=>true,
				:Lower=>1e-3
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop Correction Factor for Unequal Baffle Spacing",
				:Hidden=>true,
				:Lower=>1e-3
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Mass Velocity in Window Zone",
				:Hidden=>true,
				:Unit=>"kg/m^2/s"
			)),
			constant(Dict{Symbol,Any}(
				:Brief=>"Shell Side Ji Factor",
				:Symbol=>"J_i",
				:Hidden=>true,
				:Default=>0.05
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Shell Side Jr Factor",
				:Symbol=>"J_r",
				:Hidden=>true,
				:Lower=>10e-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Shell Side Jl Factor",
				:Symbol=>"J_l",
				:Hidden=>true,
				:Lower=>10e-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Shell Side Jb Factor",
				:Symbol=>"J_b",
				:Hidden=>true,
				:Lower=>10e-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Shell Side Jc Factor",
				:Symbol=>"J_c",
				:Hidden=>true,
				:Lower=>10e-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Shell Side Js Factor",
				:Symbol=>"J_s",
				:Hidden=>true,
				:Lower=>10e-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Shell Side Jtotal Factor",
				:Symbol=>"J_{total}",
				:Hidden=>true,
				:Lower=>10e-6
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Shell Side Cross Flow Area",
				:Symbol=>"S_m",
				:Hidden=>true,
				:Default=>0.05,
				:Lower=>10e-6
			)),
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Components"
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight"
			)),(NComp)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Tube Side Flow Regime ",
				:Hidden=>true,
				:Valid=>["laminar","transition","turbulent"],
				:Default=>"laminar"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Shell Side Flow Regime ",
				:Hidden=>true,
				:Valid=>["deep laminar","laminar","turbulent"],
				:Default=>"deep laminar"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Shell Side Flow Regime Range for Correction Factor",
				:Hidden=>true,
				:Valid=>["range1","range2","range3", "range4","range5"],
				:Default=>"range1"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Flag for Fluid Alocation ",
				:Hidden=>true,
				:Valid=>["shell","tubes"],
				:Default=>"shell"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Tube Heat Transfer Correlation in Laminar Flow",
				:Valid=>["Hausen","Schlunder"],
				:Default=>"Hausen"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Tube Heat Transfer Correlation in Transition Flow",
				:Valid=>["Gnielinski","ESDU"],
				:Default=>"Gnielinski"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Tube Heat Transfer Correlation in Turbulent Flow",
				:Valid=>["Petukhov","SiederTate"],
				:Default=>"Petukhov"
			)),
			constant(Dict{Symbol,Any}(
				:Brief=>"Pi Number",
				:Hidden=>true,
				:Default=>3.14159265,
				:Symbol=>"\\pi"
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Shell Outlet Nozzle Area",
				:Hidden=>true,
				:Lower=>1E-6 ,
				:Symbol=>"A_{nozzle\\_out }^{shell}"
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Shell Inlet Nozzle Area",
				:Hidden=>true,
				:Lower=>1E-6 ,
				:Symbol=>"A_{nozzle\\_in }^{shell}"
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Shell Outlet Escape Area Under Nozzle",
				:Hidden=>true,
				:Lower=>1E-6 ,
				:Symbol=>"Aescape_{nozzle\\_out }^{shell}"
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Shell Inlet Escape Area Under Nozzle",
				:Hidden=>true,
				:Lower=>1E-6 ,
				:Symbol=>"Aescape_{nozzle\\_in }^{shell}"
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Tube Outlet Nozzle Area",
				:Hidden=>true,
				:Lower=>1E-6 ,
				:Symbol=>"A_{nozzle\\_out }^{tube}"
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Tube Inlet Nozzle Area",
				:Hidden=>true,
				:Lower=>1E-6 ,
				:Symbol=>"A_{nozzle\\_in }^{tube}"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Tube Inlet Nozzle Pressure Loss Coeff",
				:Hidden=>true,
				:Default=>1.1,
				:Symbol=>"K_{in }^{tube}"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Tube Outlet Nozzle Pressure Loss Coeff",
				:Hidden=>true,
				:Default=>0.7,
				:Symbol=>"K_{out }^{tube}"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Baffle cut angle in degrees",
				:Symbol=>"\\theta _{ds}",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Baffle cut angle relative to the centerline in degrees",
				:Symbol=>"\\theta _{ctl}",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Fraction of number of tubes in baffle window",
				:Symbol=>"F _{tw}",
				:Hidden=>true
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Shell to baffle leakage area",
				:Symbol=>"S _{cd}",
				:Hidden=>true
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Tube to baffle hole leakage area",
				:Symbol=>"S _{td}",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Ratio of the shell to baffle leakage area",
				:Symbol=>"R_s",
				:Hidden=>true
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Hydraulic diameter of the baffle window",
				:Symbol=>"D _w",
				:Hidden=>true
			)),
			[
				:(Shell.Properties.Average.T = 0.5*InletShell.T + 0.5*OutletShell.T),
				:(Tubes.Properties.Average.T = 0.5*InletTube.T + 0.5*OutletTube.T),
				:(Shell.Properties.Average.P = 0.5*InletShell.P+0.5*OutletShell.P),
				:(Tubes.Properties.Average.P = 0.5*InletTube.P+0.5*OutletTube.P),
				:(Shell.Properties.Average.Mw = sum(M*InletShell.z)),
				:(Tubes.Properties.Average.Mw = sum(M*OutletTube.z)),
				:(Tubes.Properties.Average.Cp = PP.LiquidCp(Tubes.Properties.Average.T,Tubes.Properties.Average.P,OutletTube.z)),
				:(Tubes.Properties.Average.rho = PP.LiquidDensity(Tubes.Properties.Average.T,Tubes.Properties.Average.P,OutletTube.z)),
				:(Tubes.Properties.Inlet.rho = PP.LiquidDensity(OutletTube.T,OutletTube.P,OutletTube.z)),
				:(Tubes.Properties.Outlet.rho = PP.LiquidDensity(OutletTube.T,OutletTube.P,OutletTube.z)),
				:(Tubes.Properties.Average.Mu = PP.LiquidViscosity(Tubes.Properties.Average.T,Tubes.Properties.Average.P,OutletTube.z)),
				:(Tubes.Properties.Average.K = PP.LiquidThermalConductivity(Tubes.Properties.Average.T,Tubes.Properties.Average.P,OutletTube.z)),
				:(Tubes.Properties.Wall.Mu = PP.LiquidViscosity(Tubes.Properties.Wall.Twall,Tubes.Properties.Average.P,OutletTube.z)),
				:(Tubes.Properties.Average.Cp = PP.VapourCp(Tubes.Properties.Average.T,Tubes.Properties.Average.P,OutletTube.z)),
				:(Tubes.Properties.Average.rho = PP.VapourDensity(Tubes.Properties.Average.T,Tubes.Properties.Average.P,OutletTube.z)),
				:(Tubes.Properties.Inlet.rho = PP.VapourDensity(OutletTube.T,OutletTube.P,OutletTube.z)),
				:(Tubes.Properties.Outlet.rho = PP.VapourDensity(OutletTube.T,OutletTube.P,OutletTube.z)),
				:(Tubes.Properties.Average.Mu = PP.VapourViscosity(Tubes.Properties.Average.T,Tubes.Properties.Average.P,OutletTube.z)),
				:(Tubes.Properties.Average.K = PP.VapourThermalConductivity(Tubes.Properties.Average.T,Tubes.Properties.Average.P,OutletTube.z)),
				:(Tubes.Properties.Wall.Mu = PP.VapourViscosity(Tubes.Properties.Wall.Twall,Tubes.Properties.Average.P,OutletTube.z)),
				:(Shell.Properties.Average.Cp = PP.LiquidCp(Shell.Properties.Average.T,Shell.Properties.Average.P,InletShell.z)),
				:(Shell.Properties.Average.rho = PP.LiquidDensity(Shell.Properties.Average.T,Shell.Properties.Average.P,InletShell.z)),
				:(Shell.Properties.Inlet.rho = PP.LiquidDensity(InletShell.T,InletShell.P,InletShell.z)),
				:(Shell.Properties.Outlet.rho = PP.LiquidDensity(OutletShell.T,OutletShell.P,OutletShell.z)),
				:(Shell.Properties.Average.Mu = PP.LiquidViscosity(Shell.Properties.Average.T,Shell.Properties.Average.P,InletShell.z)),
				:(Shell.Properties.Average.K = PP.LiquidThermalConductivity(Shell.Properties.Average.T,Shell.Properties.Average.P,InletShell.z)),
				:(Shell.Properties.Wall.Mu = PP.LiquidViscosity(Shell.Properties.Wall.Twall,Shell.Properties.Average.P,InletShell.z)),
				:(Shell.Properties.Average.Cp = PP.VapourCp(Shell.Properties.Average.T,Shell.Properties.Average.P,InletShell.z)),
				:(Shell.Properties.Average.rho = PP.VapourDensity(Shell.Properties.Average.T,Shell.Properties.Average.P,InletShell.z)),
				:(Shell.Properties.Inlet.rho = PP.VapourDensity(InletShell.T,InletShell.P,InletShell.z)),
				:(Shell.Properties.Outlet.rho = PP.VapourDensity(OutletShell.T,OutletShell.P,OutletShell.z)),
				:(Shell.Properties.Average.Mu = PP.VapourViscosity(Shell.Properties.Average.T,Shell.Properties.Average.P,InletShell.z)),
				:(Shell.Properties.Average.K = PP.VapourThermalConductivity(Shell.Properties.Average.T,Shell.Properties.Average.P,InletShell.z)),
				:(Shell.Properties.Wall.Mu = PP.VapourViscosity(Shell.Properties.Wall.Twall,Shell.Properties.Average.P,InletShell.z)),
				:(Details.Q = InletShell.F*(InletShell.h-OutletShell.h)),
				:(Details.Q =-InletTube.F*(InletTube.h-OutletTube.h)),
				:(Details.Q = InletTube.F*(InletTube.h-OutletTube.h)),
				:(Details.Q =-InletShell.F*(InletShell.h-OutletShell.h)),
				:(Tubes.Properties.Inlet.Fw = sum(M*InletTube.z)*InletTube.F),
				:(Tubes.Properties.Outlet.Fw = sum(M*OutletTube.z)*OutletTube.F),
				:(Shell.Properties.Inlet.Fw = sum(M*InletShell.z)*InletShell.F),
				:(Shell.Properties.Outlet.Fw = sum(M*OutletShell.z)*OutletShell.F),
				:(OutletShell.F = InletShell.F),
				:(OutletTube.F = InletTube.F),
				:(OutletShell.z=InletShell.z),
				:(OutletTube.z=InletTube.z),
				:(Jc = 0.55+0.72*(1-2*Ftw)),
				:(Jl = 0.44*(1-Rs)+(1-0.44*(1-Rs))*exp(-2.2*(Scd + Std)/Sm)),
				:(Jtotal = Jc*Jl*Jb*Jr*Js),
				:(mw = Shell.Properties.Inlet.Fw/sqrt(abs(Sm*abs((Pi*Shell.ShellID*Shell.ShellID*((Ods/360)-sin(Ods*Pi/180)/(2*Pi))/4)-(Tubes.NumberOfTubes*Pi*Tubes.TubeOD*Tubes.TubeOD*Ftw/4))))),
				:(Tubes.HeatTransfer.fi = 16/Tubes.HeatTransfer.Re),
				:(Tubes.PressureDrop.FricFactor = 16/Tubes.HeatTransfer.Re),
				:(Tubes.HeatTransfer.Nu = 3.665 + ((0.19*((Tubes.TubeID/Tubes.TubeLength)*Tubes.HeatTransfer.Re*Tubes.HeatTransfer.PR)^0.8)/(1+0.117*((Tubes.TubeID/Tubes.TubeLength)*Tubes.HeatTransfer.Re*Tubes.HeatTransfer.PR)^0.467))),
				:(Tubes.HeatTransfer.Nu = (49.027896+4.173281*Tubes.HeatTransfer.Re*Tubes.HeatTransfer.PR*(Tubes.TubeID/Tubes.TubeLength))^(1/3)),
				:(Tubes.HeatTransfer.fi = 1/(0.79*ln(Tubes.HeatTransfer.Re)-1.64)^2),
				:(Tubes.PressureDrop.FricFactor = 0.0122),
				:(Tubes.HeatTransfer.Nu*(1+(12.7*sqrt(0.125*Tubes.HeatTransfer.fi)*((Tubes.HeatTransfer.PR)^(2/3) -1))) = 0.125*Tubes.HeatTransfer.fi*(Tubes.HeatTransfer.Re-1000)*Tubes.HeatTransfer.PR),
				:(Tubes.HeatTransfer.Nu =1),
				:(Tubes.HeatTransfer.fi = 1/(1.82*log(Tubes.HeatTransfer.Re)-1.64)^2),
				:(Tubes.PressureDrop.FricFactor = 0.0035 + 0.264*Tubes.HeatTransfer.Re^(-0.42)),
				:(Tubes.HeatTransfer.Nu*(1.07+(12.7*sqrt(0.125*Tubes.HeatTransfer.fi)*((Tubes.HeatTransfer.PR)^(2/3) -1))) = 0.125*Tubes.HeatTransfer.fi*Tubes.HeatTransfer.Re*Tubes.HeatTransfer.PR),
				:(Tubes.HeatTransfer.Nu = 0.027*(Tubes.HeatTransfer.PR)^(1/3)*(Tubes.HeatTransfer.Re)^(4/5)),
				:(Sm= Baffles.Central_Spacing*(Clearances.BundleToShell+((Shell.ShellID-Clearances.BundleToShell-Tubes.TubeOD)/Tubes.TubePitch)*(Tubes.TubePitch-Tubes.TubeOD))),
				:(Nc = Shell.ShellID*(1-0.02*Baffles.BaffleCut)/(0.866*Tubes.TubePitch)),
				:(Ncw = 0.8*(Shell.ShellID*0.01*Baffles.BaffleCut-(Clearances.BundleToShell + Tubes.TubeOD)*0.5)/(0.866*Tubes.TubePitch)),
				:(a = 1.45/(1+0.14*Shell.HeatTransfer.Re^0.519)),
				:(b=7/(1+0.14*Shell.HeatTransfer.Re^0.5)),
				:(Rss = Clearances.SealStrip/(Shell.ShellID*(1-0.02*Baffles.BaffleCut)/(0.866*Tubes.TubePitch))),
				:(Shell.PressureDrop.Ideal= 2*Shell.PressureDrop.FricFactor*(Shell.ShellID*(1-0.02*Baffles.BaffleCut)/(0.866*Tubes.TubePitch))*(Shell.Properties.Inlet.Fw/Sm)^2/(Shell.Properties.Average.rho*Shell.HeatTransfer.Phi)),
				:(Shell.PressureDrop.EndZones = Shell.PressureDrop.Ideal*(1+ (Ncw/(Shell.ShellID*(1-0.02*Baffles.BaffleCut)/(0.866*Tubes.TubePitch))))*Rb*Rspd),
				:(Ji =1.40*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^0.667),
				:(Shell.PressureDrop.FricFactor=48*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-1)),
				:(Ji =1.36*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^(-0.657)),
				:(Shell.PressureDrop.FricFactor=45.10*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-0.973)),
				:(Ji =0.593*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^(-0.477)),
				:(Shell.PressureDrop.FricFactor=4.570*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-0.476)),
				:(Ji =0.321*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^(-0.388)),
				:(Shell.PressureDrop.FricFactor=0.486*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-0.152)),
				:(Ji =0.321*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^(-0.388)),
				:(Shell.PressureDrop.FricFactor=0.372*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-0.123)),
				:(Sm= Baffles.Central_Spacing*(Clearances.BundleToShell+((Shell.ShellID-Clearances.BundleToShell-Tubes.TubeOD)/(0.707*Tubes.TubePitch))*(Tubes.TubePitch-Tubes.TubeOD))),
				:(Nc = Shell.ShellID*(1-0.02*Baffles.BaffleCut)/(0.707*Tubes.TubePitch)),
				:(Ncw = 0.8*(Shell.ShellID*0.01*Baffles.BaffleCut-(Clearances.BundleToShell + Tubes.TubeOD)*0.5)/(0.707*Tubes.TubePitch)),
				:(a = 1.930/(1+0.14*Shell.HeatTransfer.Re^0.500)),
				:(b=6.59/(1+0.14*Shell.HeatTransfer.Re^0.52)),
				:(Rss = Clearances.SealStrip/(Shell.ShellID*(1-0.02*Baffles.BaffleCut)/(0.707*Tubes.TubePitch))),
				:(Shell.PressureDrop.Ideal= 2*Shell.PressureDrop.FricFactor*(Shell.ShellID*(1-0.02*Baffles.BaffleCut)/(0.707*Tubes.TubePitch))*(Shell.Properties.Inlet.Fw/Sm)^2/(Shell.Properties.Average.rho*Shell.HeatTransfer.Phi)),
				:(Shell.PressureDrop.EndZones = Shell.PressureDrop.Ideal*(1+ (Ncw/(Shell.ShellID*(1-0.02*Baffles.BaffleCut)/(0.707*Tubes.TubePitch))))*Rb*Rspd),
				:(Ji =1.550*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^0.667),
				:(Shell.PressureDrop.FricFactor=32*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-1)),
				:(Ji =0.498*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^0.656),
				:(Shell.PressureDrop.FricFactor=26.20*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-0.913)),
				:(Ji =0.730*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^0.500),
				:(Shell.PressureDrop.FricFactor=3.50*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-0.476)),
				:(Ji =0.370*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^(-0.396)),
				:(Shell.PressureDrop.FricFactor=0.333*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-0.136)),
				:(Ji =0.370*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^(-0.396)),
				:(Shell.PressureDrop.FricFactor=0.303*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-0.126)),
				:(Sm= Baffles.Central_Spacing*(Clearances.BundleToShell+((Shell.ShellID-Clearances.BundleToShell-Tubes.TubeOD)/Tubes.TubePitch)*(Tubes.TubePitch-Tubes.TubeOD))),
				:(Nc = Shell.ShellID*(1-0.02*Baffles.BaffleCut)/Tubes.TubePitch),
				:(Ncw = 0.8*(Shell.ShellID*0.01*Baffles.BaffleCut-(Clearances.BundleToShell + Tubes.TubeOD)*0.5)/Tubes.TubePitch),
				:(a = 1.187/(1+0.14*Shell.HeatTransfer.Re^0.370)),
				:(b=6.30/(1+0.14*Shell.HeatTransfer.Re^0.38)),
				:(Rss = Clearances.SealStrip/(Shell.ShellID*(1-0.02*Baffles.BaffleCut)/Tubes.TubePitch)),
				:(Shell.PressureDrop.Ideal= 2*Shell.PressureDrop.FricFactor*(Shell.ShellID*(1-0.02*Baffles.BaffleCut)/Tubes.TubePitch)*(Shell.Properties.Inlet.Fw/Sm)^2/(Shell.Properties.Average.rho*Shell.HeatTransfer.Phi)),
				:(Shell.PressureDrop.EndZones = Shell.PressureDrop.Ideal*(1+ (Ncw/(Shell.ShellID*(1-0.02*Baffles.BaffleCut)/Tubes.TubePitch)))*Rb*Rspd),
				:(Ji =0.970*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^(-0.667)),
				:(Shell.PressureDrop.FricFactor=35*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-1)),
				:(Ji =0.900*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^(-0.631)),
				:(Shell.PressureDrop.FricFactor=32.10*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-0.963)),
				:(Ji =0.408*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^(-0.460)),
				:(Shell.PressureDrop.FricFactor=6.090*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-0.602)),
				:(Ji =0.107*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^(-0.266)),
				:(Shell.PressureDrop.FricFactor=0.0815*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^0.022),
				:(Ji =0.370*((1.33*Tubes.TubeOD/Tubes.TubePitch)^a)*Shell.HeatTransfer.Re^(-0.395)),
				:(Shell.PressureDrop.FricFactor=0.391*((1.33*Tubes.TubeOD/Tubes.TubePitch)^b)*Shell.HeatTransfer.Re^(-0.148)),
				:(Jr = (10/((Nc +Ncw)*(Baffles.NumberOfBaffles+1)))^0.18),
				:(Js = (Baffles.NumberOfBaffles-1+(Baffles.Inlet_Spacing/Baffles.Central_Spacing)^0.7 + (Baffles.Outlet_Spacing/Baffles.Central_Spacing)^0.7)/(Baffles.NumberOfBaffles-1+(Baffles.Inlet_Spacing/Baffles.Central_Spacing) + (Baffles.Outlet_Spacing/Baffles.Central_Spacing))),
				:(Jb = exp(-1.35*( Clearances.BundleToShell+ Tubes.TubeOD)*Baffles.Central_Spacing/Sm*(1-(2*(Clearances.SealStrip/Nc)^(1/3))))),
				:(Rb = exp(-4.7*((Clearances.BundleToShell + Tubes.TubeOD)*Baffles.Central_Spacing/Sm)*(1-(2*Rss)^(1/3)))),
				:(Rspd = (Baffles.Central_Spacing/Baffles.Outlet_Spacing) + (Baffles.Central_Spacing/Baffles.Inlet_Spacing)),
				:(Shell.PressureDrop.Window = Baffles.NumberOfBaffles*((26/Shell.Properties.Average.rho)*mw*Shell.Properties.Average.Mu*(Ncw/(Tubes.TubePitch-Tubes.TubeOD)+ Baffles.Central_Spacing/(Dw*Dw))+ 0.5*mw*mw/Shell.Properties.Average.rho)*exp(-1.33*(1+Rs)*((Scd + Std)/Sm)^(-0.15*(1+Rs) + 0.8))),
				:(Jr = (10/((Nc +Ncw)*(Baffles.NumberOfBaffles+1)))^0.18 + (0.25-0.0125*Shell.HeatTransfer.Re)*((10/((Nc +Ncw)*(Baffles.NumberOfBaffles+1)))^0.18 - 1)),
				:(Js = (Baffles.NumberOfBaffles-1+(Baffles.Inlet_Spacing/Baffles.Central_Spacing)^0.7 + (Baffles.Outlet_Spacing/Baffles.Central_Spacing)^0.7)/(Baffles.NumberOfBaffles-1+(Baffles.Inlet_Spacing/Baffles.Central_Spacing) + (Baffles.Outlet_Spacing/Baffles.Central_Spacing))),
				:(Jb = exp(-1.35*( Clearances.BundleToShell+ Tubes.TubeOD)*Baffles.Central_Spacing/Sm*(1-(2*(Clearances.SealStrip/Nc)^(1/3))))),
				:(Rb = exp(-4.7*((Clearances.BundleToShell+ Tubes.TubeOD)*Baffles.Central_Spacing/Sm)*(1-(2*Rss)^(1/3)))),
				:(Rspd = (Baffles.Central_Spacing/Baffles.Outlet_Spacing) + (Baffles.Central_Spacing/Baffles.Inlet_Spacing)),
				:(Shell.PressureDrop.Window = Baffles.NumberOfBaffles*((26/Shell.Properties.Average.rho)*mw*Shell.Properties.Average.Mu*(Ncw/(Tubes.TubePitch-Tubes.TubeOD)+ Baffles.Central_Spacing/(Dw*Dw))+ 0.5*mw*mw/Shell.Properties.Average.rho)*exp(-1.33*(1+Rs)*((Scd + Std)/Sm)^(-0.15*(1+Rs) + 0.8))),
				:(Jr = 1),
				:(Js = (Baffles.NumberOfBaffles-1+(Baffles.Inlet_Spacing/Baffles.Central_Spacing)^0.4 + (Baffles.Outlet_Spacing/Baffles.Central_Spacing)^0.4)/(Baffles.NumberOfBaffles-1+(Baffles.Inlet_Spacing/Baffles.Central_Spacing) + (Baffles.Outlet_Spacing/Baffles.Central_Spacing))),
				:(Jb = exp(-1.25*( Clearances.BundleToShell+ Tubes.TubeOD)*Baffles.Central_Spacing/Sm*(1-(2*(Clearances.SealStrip/Nc)^(1/3))))),
				:(Rb = exp(-3.7*((Clearances.BundleToShell + Tubes.TubeOD)*Baffles.Central_Spacing/Sm)*(1-(2*Rss)^(1/3)))),
				:(Rspd = (Baffles.Central_Spacing/Baffles.Outlet_Spacing)^1.8 + (Baffles.Central_Spacing/Baffles.Inlet_Spacing)^1.8),
				:(Shell.PressureDrop.Window = Baffles.NumberOfBaffles*((2+0.6*Ncw)*0.5*mw*mw/Shell.Properties.Average.rho)*exp(-1.33*(1+Rs)*((Scd + Std)/Sm)^(-0.15*(1+Rs) + 0.8))),
				:(Shell.PressureDrop.CrossFlow= Shell.PressureDrop.Ideal*Rb*(Baffles.NumberOfBaffles-1)*exp(-1.33*(1+Rs)*((Scd + Std)/Sm)^(-0.15*(1+Rs) + 0.8))),
				:(Shell.HeatTransfer.Phi = (Shell.Properties.Average.Mu/Shell.Properties.Wall.Mu)^0.14),
				:(Tubes.HeatTransfer.Phi = (Tubes.Properties.Average.Mu/Tubes.Properties.Wall.Mu)^0.14),
				:(Shell.PressureDrop.RVsquare_in = Shell.Properties.Inlet.rho*(Shell.PressureDrop.Vnozzle_in)^2),
				:(Shell.PressureDrop.RVsquare_out = Shell.Properties.Outlet.rho*(Shell.PressureDrop.Vnozzle_out)^2),
				:(Tubes.PressureDrop.TubeFriction = 2*Tubes.PressureDrop.FricFactor*Tubes.TubeLength*Tubes.Properties.Average.rho*(Tubes.HeatTransfer.Vtube^2)*Tubes.Tubepasses/(Tubes.TubeID*Tubes.HeatTransfer.Phi)),
				:(Tubes.PressureDrop.InletNozzle = 0.5*Kinlet_Tube*Tubes.Properties.Inlet.rho*Tubes.PressureDrop.Vnozzle_in^2),
				:(Tubes.PressureDrop.Vnozzle_in = Tubes.Properties.Inlet.Fw/(Tubes.Properties.Inlet.rho*Ainozzle_Tube)),
				:(Tubes.PressureDrop.OutletNozzle = 0.5*Koutlet_Tube*Tubes.Properties.Outlet.rho*Tubes.PressureDrop.Vnozzle_out^2),
				:(Tubes.PressureDrop.Vnozzle_out = Tubes.Properties.Inlet.Fw/(Tubes.Properties.Outlet.rho*Aonozzle_Tube)),
				:(Shell.PressureDrop.InletNozzle = (0.5*Shell.Properties.Inlet.Fw^2/Shell.Properties.Inlet.rho)*((1/Ainozzle_Shell^2)+(1/Aeinozzle_Shell^2))),
				:(Shell.PressureDrop.Vnozzle_in = Shell.Properties.Inlet.Fw/(Shell.Properties.Inlet.rho*Ainozzle_Shell)),
				:(Shell.PressureDrop.OutletNozzle = (0.5*Shell.Properties.Outlet.Fw^2/Shell.Properties.Outlet.rho)*((1/Ainozzle_Shell^2)+(1/Aeinozzle_Shell^2))),
				:(Shell.PressureDrop.Vnozzle_out = Shell.Properties.Outlet.Fw/(Shell.Properties.Outlet.rho*Aonozzle_Shell)),
				:(OutletShell.P = InletShell.P - Shell.PressureDrop.Total),
				:(OutletTube.P = InletTube.P - Tubes.PressureDrop.Total),
				:(Shell.Properties.Wall.Twall = (Shell.Properties.Average.T+Tubes.Properties.Average.T)/2),
				:(Tubes.Properties.Wall.Twall = (Shell.Properties.Average.T+Tubes.Properties.Average.T)/2),
				:(Tubes.HeatTransfer.Vtube = Tubes.Properties.Inlet.Fw*Tubes.Tubepasses/((Pi*Tubes.TubeID*Tubes.TubeID/4)*Tubes.Properties.Average.rho*Tubes.NumberOfTubes)),
				:(Tubes.HeatTransfer.Re = (Tubes.Properties.Average.rho*Tubes.HeatTransfer.Vtube*Tubes.TubeID)/Tubes.Properties.Average.Mu),
				:(Tubes.HeatTransfer.PR = ((Tubes.Properties.Average.Cp/Tubes.Properties.Average.Mw)*Tubes.Properties.Average.Mu)/Tubes.Properties.Average.K),
				:(Tubes.HeatTransfer.htube= (Tubes.HeatTransfer.Nu*Tubes.Properties.Average.K/Tubes.TubeID)*Tubes.HeatTransfer.Phi),
				:(Shell.HeatTransfer.PR = ((Shell.Properties.Average.Cp/Shell.Properties.Average.Mw)*Shell.Properties.Average.Mu)/Shell.Properties.Average.K),
				:(Details.Ud=1/(Tubes.TubeOD/(Tubes.HeatTransfer.htube*Tubes.TubeID)+Shell.Fouling+Tubes.Fouling*(Tubes.TubeOD/Tubes.TubeID)+(Tubes.TubeOD*ln(Tubes.TubeOD/Tubes.TubeID)/(2*Tubes.Kwall))+(1/(Shell.HeatTransfer.hshell)))),
				:(Details.Uc=1/(Tubes.TubeOD/(Tubes.HeatTransfer.htube*Tubes.TubeID)+(Tubes.TubeOD*ln(Tubes.TubeOD/Tubes.TubeID)/(2*Tubes.Kwall))+(1/(Shell.HeatTransfer.hshell)))),
				:(Details.A=Pi*Tubes.TubeOD*Tubes.NumberOfTubes*Tubes.TubeLength),
				:(Tubes.TubeLength = Baffles.Inlet_Spacing+Baffles.Outlet_Spacing+Baffles.Central_Spacing*(Baffles.NumberOfBaffles-1)),
				:(Shell.HeatTransfer.Re = (Tubes.TubeOD*Shell.Properties.Inlet.Fw/Sm)/Shell.Properties.Average.Mu),
				:(Shell.HeatTransfer.hshell = Ji*(Shell.Properties.Average.Cp/Shell.Properties.Average.Mw)*(Shell.Properties.Inlet.Fw/Sm)*(Shell.HeatTransfer.PR^(-2/3))*Jtotal*Shell.HeatTransfer.Phi),
			],
			[
				"Shell Stream Average Temperature","Tube Stream Average Temperature","Shell Stream Average Pressure","Tube Stream Average Pressure","Shell Stream Average Molecular Weight","Tube Stream Average Molecular Weight","Tube Stream Average Heat Capacity","Tube Stream Average Mass Density","Tube Stream Inlet Mass Density","Tube Stream Outlet Mass Density","TubeStream Average Viscosity","Tube Stream Average Conductivity","Tube Stream Viscosity at Wall Temperature","Tube Stream Average Heat Capacity","Tube Stream Average Mass Density","Tube Stream Inlet Mass Density","Tube Stream Outlet Mass Density","Tube Stream Average Viscosity ","Tube Stream Average Conductivity ","Tube Stream Viscosity at Wall Temperature","Shell Stream Average Heat Capacity","Shell Stream Average Mass Density","ShellStream Inlet Mass Density","Shell Stream Outlet Mass Density","Shell Stream Average Viscosity","Shell Stream Average Conductivity","ShellStream Viscosity  at Wall Temperature","Shell Stream Average Heat Capacity","Shell Stream Average Mass Density","Shell Stream Inlet Mass Density","Shell Stream Outlet Mass Density","Shell Stream Average Viscosity","Shell Stream Average Conductivity","Shell Stream Viscosity at Wall Temperature","Energy Balance Hot Stream","Energy Balance Cold Stream","Energy Balance Hot Stream","Energy Balance Cold Stream","Flow Mass Inlet Tube Stream","Flow Mass Outlet Tube Stream","Flow Mass Inlet Shell Stream","Flow Mass Outlet Shell Stream","Molar Balance Shell Stream","Molar Balance Tube Stream","Shell	Stream Molar Fraction Constraint","Tube Stream Molar Fraction Constraint","Jc Factor","Jl Factor","Total J Factor","Mass Velocity in Window Zone","Friction Factor for heat Transfer: Not Necessary in Laminar Correlation - Use any one equation that you want","Friction Factor for Pressure Drop in Laminar Flow","Nusselt Number in Laminar Flow - Hausen Equation","Nusselt Number in Laminar Flow - Schlunder Equation","Friction Factor for heat Transfer : for use in Gnielinski Equation","Friction Factor for Pressure Drop in Transition Flow","Nusselt Number in Transition Flow - Gnielinski Equation","Nusselt Number in Transition Flow - ESDU Equation","Friction Factor for heat Transfer : for use in Petukhov Equation","Friction Factor for Pressure Drop in Turbulent Flow","Nusselt Number in Turbulent Flow - Petukhov Equation","Nusselt Number in Transition Flow - Sieder Tate Equation","Shell Side Cross Flow Area","Number of Tube rows Crossed in one Crossflow Section","Number of Effective Crossflow rows in Each Window","Variable for calculating Ji heat transfer correction Factor","Variable for calculating Shell Side Pressure Drop Friction Factor","Correction Factor for Pressure Drop","Ideal Shell Side Pressure Drop","Shell Pressure End Zones","Ji Factor","Shell Side Pressure Drop Friction Factor","Ji Factor","Shell Side Pressure Drop Friction Factor","Ji Factor","Shell Side Pressure Drop Friction Factor","Ji Factor","Shell Side Pressure Drop Friction Factor","Ji Factor","Shell Side Pressure Drop Friction Factor","Shell Side Cross Flow Area","Number of Tube rows Crossed in one Crossflow Section","Number of Effective Crossflow rows in Each Window","Variable for calculating Ji heat transfer correction Factor","Variable for calculating Shell Side Pressure Drop Friction Factor","Correction Factor for Pressure Drop","Ideal Shell Side Pressure Drop","Shell Pressure End Zones","Ji Factor","Shell Side Pressure Drop Friction Factor","Ji Factor","Shell Side Pressure Drop Friction Factor","Ji Factor","Shell Side Pressure Drop Friction Factor","Ji Factor","Shell Side Pressure Drop Friction Factor","Ji Factor","Shell Side Pressure Drop Friction Factor","Shell Side Cross Flow Area","Number of Tube rows Crossed in one Crossflow Section","Number of Effective Crossflow rows in Each Window","Variable for calculating Ji heat transfer correction Factor","Variable for calculating Shell Side Pressure Drop Friction Factor","Correction Factor for Pressure Drop","Ideal Shell Side Pressure Drop","Shell Pressure End Zones","Ji Factor","Shell Side Pressure Drop Friction Factor","Ji Factor","Shell Side Pressure Drop Friction Factor","Ji Factor","Shell Side Pressure Drop Friction Factor","Ji Factor","Shell Side Pressure Drop Friction Factor","Ji Factor","Shell Side Pressure Drop Friction Factor","Jr Factor","Js Factor","Jb Factor","ByPass Correction Factor for Pressure Drop","Pressure Drop Correction Factor for Unequal Baffle Spacing","Shell Pressure Drop Baffle Window","Jr Factor","Js Factor","Jb Factor","ByPass Correction Factor for Pressure Drop","Pressure Drop Correction Factor for Unequal Baffle Spacing","Shell Pressure Drop Baffle Window","Jr Factor","Js Factor","Jb Factor","ByPass Correction Factor for Pressure Drop","Pressure Drop Correction Factor for Unequal Baffle Spacing","Shell Pressure Drop Baffle Window","Shell Pressure Drop Cross Flow","Shell Side Phi correction","Tube Side Phi correction","Shell Side inlet Nozzle rho-V^2","Shell Side Outlet Nozzle rho-V^2","Tube Side Pressure Drop","Pressure Drop Tube Side Inlet Nozzle","Velocity Tube Side Inlet Nozzle","Pressure Drop Tube Side Outlet Nozzle","Velocity Tube Side Outlet Nozzle","Shell Pressure Drop Inlet Nozzle","Velocity Shell Side Inlet Nozzle","Shell Pressure Drop Outlet Nozzle","Velocity Shell Side Outlet Nozzle","Pressure Drop Shell Stream","Pressure Drop Tube Stream","Shell Wall Temperature","Tube Wall Temperature","Tube Side Velocity","Tube Side Reynolds Number","Tube Side Prandtl Number","Tube Side Film Coefficient","Shell Side Prandtl Number","Overall Heat Transfer Coefficient Dirty","Overall Heat Transfer Coefficient Clean","Exchange Surface Area","Baffle Spacing Constraint","Shell Side Reynolds Number","Shell Heat Transfer Coefficient",
			],
			[:HotSide,:ShellType,:Pattern,:PP,:NComp,:M,:TubeFlowRegime,:ShellFlowRegime,:ShellRange,:Side,:LaminarCorrelation,:TransitionCorrelation,:TurbulentCorrelation,:Pi,:Aonozzle_Shell,:Ainozzle_Shell,:Aeonozzle_Shell,:Aeinozzle_Shell,:Aonozzle_Tube,:Ainozzle_Tube,:Kinlet_Tube,:Koutlet_Tube,:Ods,:Octl,:Ftw,:Scd,:Std,:Rs,:Dw,],
			[:Tubes,:Shell,:Baffles,:Clearances,:InletTube,:OutletTube,:InletShell,:OutletShell,:Details,:Nc,:Ncw,:a,:b,:Rb,:Rss,:Rspd,:mw,:Ji,:Jr,:Jl,:Jb,:Jc,:Js,:Jtotal,:Sm,]
		)
	end
	HotSide::DanaSwitcher
	ShellType::DanaSwitcher
	Pattern::DanaSwitcher
	Tubes::Tube_Side_Main
	Shell::Shell_Side_Main
	Baffles::Baffles_Main
	Clearances::Clearances_Main
	InletTube::stream
	OutletTube::streamPH
	InletShell::stream
	OutletShell::streamPH
	Details::Details_Main
	Nc::DanaReal
	Ncw::DanaReal
	a::DanaReal
	b::DanaReal
	Rb::DanaReal
	Rss::DanaReal
	Rspd::DanaReal
	mw::DanaReal
	Ji::constant
	Jr::positive
	Jl::positive
	Jb::positive
	Jc::positive
	Js::positive
	Jtotal::positive
	Sm::area
	PP::DanaPlugin
	NComp::DanaInteger
	M::Array{molweight}
	TubeFlowRegime::DanaSwitcher
	ShellFlowRegime::DanaSwitcher
	ShellRange::DanaSwitcher
	Side::DanaSwitcher
	LaminarCorrelation::DanaSwitcher
	TransitionCorrelation::DanaSwitcher
	TurbulentCorrelation::DanaSwitcher
	Pi::constant
	Aonozzle_Shell::area
	Ainozzle_Shell::area
	Aeonozzle_Shell::area
	Aeinozzle_Shell::area
	Aonozzle_Tube::area
	Ainozzle_Tube::area
	Kinlet_Tube::positive
	Koutlet_Tube::positive
	Ods::DanaReal
	Octl::DanaReal
	Ftw::DanaReal
	Scd::area
	Std::area
	Rs::DanaReal
	Dw::length
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export ShellandTubesBasic
function set(in::ShellandTubesBasic)
	#"Molecular weight"
	M = PP.MolecularWeight()
	 #"Pi Number"
	Pi = 3.14159265
	 #"Baffle cut angle in degrees"
	Ods = (360/Pi)*acos(1-0.02*Baffles.BaffleCut)
	 #"Baffle cut angle relative to the centerline"
	Octl = (360/Pi)*acos((Shell.ShellID/(Shell.ShellID - Clearances.BundleToShell - Tubes.TubeOD))*(1-0.02*Baffles.BaffleCut))
	 #"Fraction of number of tubes in baffle window"
	Ftw = (Octl/360)-sin(Octl*Pi/180)/(2*Pi)
	 #"Shell to baffle leakage area"
	Scd = Pi*Shell.ShellID*Clearances.BaffleToShell *((360-Ods)/720)
	 #"Tube to baffle hole leakage area"
	Std = Pi*0.25*((Clearances.TubeToBaffle + Tubes.TubeOD)^2-Tubes.TubeOD*Tubes.TubeOD)*Tubes.NumberOfTubes*(1-Ftw)
	 #"comments"
	Rs = Scd/(Scd+Std)
	 #"comments"
	Dw = (4*abs((Pi*Shell.ShellID*Shell.ShellID*((Ods/360)-sin(Ods*Pi/180)/(2*Pi))/4)-(Tubes.NumberOfTubes*Pi*Tubes.TubeOD*Tubes.TubeOD*Ftw/4)))/(Pi*Tubes.TubeOD*Tubes.NumberOfTubes*Ftw+ Pi*Shell.ShellID*Ods/360)
	 #"Tube Side Inlet Nozzle Area"
	Ainozzle_Tube = 0.25*Pi*Tubes.InletNozzleID^2
	 #"Tube Side Outlet Nozzle Area"
	Aonozzle_Tube = 0.25*Pi*Tubes.OutletNozzleID^2
	 #"Tube Inlet Nozzle Pressure Loss Coeff"
	Kinlet_Tube = 1.1
	 #"Tube Outlet Nozzle Pressure Loss Coeff"
	Koutlet_Tube = 0.7
	 #"Shell Outlet Nozzle Area"
	Aonozzle_Shell = 0.25*Pi*Shell.OutletNozzleID^2
	 #"Shell Inlet Nozzle Area"
	Ainozzle_Shell = 0.25*Pi*Shell.InletNozzleID^2
	 #"Shell Outlet Escape Area Under Nozzle"
	Aeonozzle_Shell = Pi*Shell.OutletNozzleID*Clearances.Honozzle_Shell + 0.6*Aonozzle_Shell*(1-(Tubes.TubeOD/Tubes.TubePitch))
	 #"Shell Inlet Escape Area Under Nozzle"
	Aeinozzle_Shell = Pi*Shell.InletNozzleID*Clearances.Hinozzle_Shell + 0.6*Ainozzle_Shell*(1-(Tubes.TubeOD/Tubes.TubePitch))

end
function setEquationFlow(in::ShellandTubesBasic)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	if InletTube.v == 0
		addEquation(7)
		addEquation(8)
		addEquation(9)
		addEquation(10)
		addEquation(11)
		addEquation(12)
		addEquation(13)
	else
		addEquation(14)
		addEquation(15)
		addEquation(16)
		addEquation(17)
		addEquation(18)
		addEquation(19)
		addEquation(20)
	end
	if InletShell.v == 0
		addEquation(21)
		addEquation(22)
		addEquation(23)
		addEquation(24)
		addEquation(25)
		addEquation(26)
		addEquation(27)
	else
		addEquation(28)
		addEquation(29)
		addEquation(30)
		addEquation(31)
		addEquation(32)
		addEquation(33)
		addEquation(34)
	end
	let switch=Side
		if Side==InletTube.T > InletShell.T
			set(switch,"tubes")
		end
		if Side==InletTube.T < InletShell.T
			set(switch,"shell")
		end
		if switch=="shell"
			addEquation(35)
			addEquation(36)
		elseif switch=="tubes"
			addEquation(37)
			addEquation(38)
		end
	end
	addEquation(39)
	addEquation(40)
	addEquation(41)
	addEquation(42)
	addEquation(43)
	addEquation(44)
	addEquation(45)
	addEquation(46)
	addEquation(47)
	addEquation(48)
	addEquation(49)
	addEquation(50)
	let switch=TubeFlowRegime
		if TubeFlowRegime==Tubes.HeatTransfer.Re > 2300
			set(switch,"transition")
		end
		if TubeFlowRegime==Tubes.HeatTransfer.Re < 2300
			set(switch,"laminar")
		end
		if TubeFlowRegime==Tubes.HeatTransfer.Re > 10000
			set(switch,"turbulent")
		end
		if TubeFlowRegime==Tubes.HeatTransfer.Re < 10000
			set(switch,"transition")
		end
		if switch=="laminar"
			addEquation(51)
			addEquation(52)
			let switch=LaminarCorrelation
				if switch=="Hausen"
					addEquation(53)
				elseif switch=="Schlunder"
					addEquation(54)
				end
			end
		elseif switch=="transition"
			addEquation(55)
			addEquation(56)
			let switch=TransitionCorrelation
				if switch=="Gnielinski"
					addEquation(57)
				elseif switch=="ESDU"
					addEquation(58)
					#to be implemented

				end
			end
		elseif switch=="turbulent"
			addEquation(59)
			addEquation(60)
			let switch=TurbulentCorrelation
				if switch=="Petukhov"
					addEquation(61)
				elseif switch=="SiederTate"
					addEquation(62)
				end
			end
		end
	end
	let switch=Pattern
		if switch=="Triangle"
			addEquation(63)
			addEquation(64)
			addEquation(65)
			addEquation(66)
			addEquation(67)
			addEquation(68)
			addEquation(69)
			addEquation(70)
			let switch=ShellRange
				if ShellRange==Shell.HeatTransfer.Re > 10
					set(switch,"range2")
				end
				if ShellRange==Shell.HeatTransfer.Re > 100
					set(switch,"range3")
				end
				if ShellRange==Shell.HeatTransfer.Re > 1000
					set(switch,"range4")
				end
				if ShellRange==Shell.HeatTransfer.Re > 10000
					set(switch,"range5")
				end
				if ShellRange==Shell.HeatTransfer.Re < 10000
					set(switch,"range4")
				end
				if switch=="range1"
					addEquation(71)
					addEquation(72)
				elseif switch=="range2"
					addEquation(73)
					addEquation(74)
				elseif switch=="range3"
					addEquation(75)
					addEquation(76)
				elseif switch=="range4"
					addEquation(77)
					addEquation(78)
				elseif switch=="range5"
					addEquation(79)
					addEquation(80)
				end
			end
		elseif switch=="Rotated Square"
			addEquation(81)
			addEquation(82)
			addEquation(83)
			addEquation(84)
			addEquation(85)
			addEquation(86)
			addEquation(87)
			addEquation(88)
			let switch=ShellRange
				if ShellRange==Shell.HeatTransfer.Re > 10
					set(switch,"range2")
				end
				if ShellRange==Shell.HeatTransfer.Re > 100
					set(switch,"range3")
				end
				if ShellRange==Shell.HeatTransfer.Re > 1000
					set(switch,"range4")
				end
				if ShellRange==Shell.HeatTransfer.Re > 10000
					set(switch,"range5")
				end
				if ShellRange==Shell.HeatTransfer.Re < 10000
					set(switch,"range4")
				end
				if switch=="range1"
					addEquation(89)
					addEquation(90)
				elseif switch=="range2"
					addEquation(91)
					addEquation(92)
				elseif switch=="range3"
					addEquation(93)
					addEquation(94)
				elseif switch=="range4"
					addEquation(95)
					addEquation(96)
				elseif switch=="range5"
					addEquation(97)
					addEquation(98)
				end
			end
		elseif switch=="Square"
			addEquation(99)
			addEquation(100)
			addEquation(101)
			addEquation(102)
			addEquation(103)
			addEquation(104)
			addEquation(105)
			addEquation(106)
			let switch=ShellRange
				if ShellRange==Shell.HeatTransfer.Re > 10
					set(switch,"range2")
				end
				if ShellRange==Shell.HeatTransfer.Re > 100
					set(switch,"range3")
				end
				if ShellRange==Shell.HeatTransfer.Re > 1000
					set(switch,"range4")
				end
				if ShellRange==Shell.HeatTransfer.Re > 10000
					set(switch,"range5")
				end
				if ShellRange==Shell.HeatTransfer.Re < 10000
					set(switch,"range4")
				end
				if switch=="range1"
					addEquation(107)
					addEquation(108)
				elseif switch=="range2"
					addEquation(109)
					addEquation(110)
				elseif switch=="range3"
					addEquation(111)
					addEquation(112)
				elseif switch=="range4"
					addEquation(113)
					addEquation(114)
				elseif switch=="range5"
					addEquation(115)
					addEquation(116)
				end
			end
		end
	end
	let switch=ShellFlowRegime
		if ShellFlowRegime==Shell.HeatTransfer.Re > 20
			set(switch,"laminar")
		end
		if ShellFlowRegime==Shell.HeatTransfer.Re < 20
			set(switch,"deep laminar")
		end
		if ShellFlowRegime==Shell.HeatTransfer.Re > 100
			set(switch,"turbulent")
		end
		if ShellFlowRegime==Shell.HeatTransfer.Re < 100
			set(switch,"laminar")
		end
		if switch=="deep laminar"
			addEquation(117)
			addEquation(118)
			addEquation(119)
			addEquation(120)
			addEquation(121)
			addEquation(122)
		elseif switch=="laminar"
			addEquation(123)
			addEquation(124)
			addEquation(125)
			addEquation(126)
			addEquation(127)
			addEquation(128)
		elseif switch=="turbulent"
			addEquation(129)
			addEquation(130)
			addEquation(131)
			addEquation(132)
			addEquation(133)
			addEquation(134)
		end
	end
	addEquation(135)
	addEquation(136)
	addEquation(137)
	addEquation(138)
	addEquation(139)
	addEquation(140)
	addEquation(141)
	addEquation(142)
	addEquation(143)
	addEquation(144)
	addEquation(145)
	addEquation(146)
	addEquation(147)
	addEquation(148)
	addEquation(149)
	addEquation(150)
	addEquation(151)
	addEquation(152)
	addEquation(153)
	addEquation(154)
	addEquation(155)
	addEquation(156)
	addEquation(157)
	addEquation(158)
	addEquation(159)
	addEquation(160)
	addEquation(161)
	addEquation(162)
	addEquation(163)
end
function atributes(in::ShellandTubesBasic,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Basic Model for Detailed Shell and Tube Heat Exchanger."
	fields[:Info]="to be documented.

== Assumptions ==
* to be documented

== Specify ==
* to be documented

== Setting Parameters ==
* to be documented

== References ==
[1] E.A.D. Saunders, Heat Exchangers: Selection, Design and
 Construction, Longman, Harlow, 1988.

[2] Taborek, J., Shell-and-tube heat exchangers, in Heat Exchanger Design Handbook, Vol. 3
 Hemisphere Publishing Corp., New York, 1988.

[3] Bell, K. J., Mueller, A. C., Wolverine Engineering Data Book II. Wolverine Tube, Inc., <www.wlv.com>, 2001.

[4] Fakheri, A. , Alternative approach for determining log mean temperature difference correction factor
 and number of shells of shell and tube heat exchangers, Journal of Enhanced Heat Transfer, v. 10, p. 407- 420, 2003.

[5] Gnielinski, V., Forced convection in ducts, in Heat Exchanger Design Handbook, Vol. 2
 Hemisphere Publishing Corp., New York, 1988."
	drive!(fields,_)
	return fields
end
ShellandTubesBasic(_::Dict{Symbol,Any})=begin
	newModel=ShellandTubesBasic()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Test
import LinearAlgebra: norm
import PhysicalConstants.CODATA2014: Z_0

function element_parabolic_phase_shift(focal_point::Vector{Float64}, element_position::Vector{Float64})
    """
    From

    """

    return phi

end

function dipole_E_field_at_point(point::Vector{Float64}, element_position::Vector{Float64}, phase_shift::Float64, current_amplitude::Float64,
                                lambda::Float64, dipole_length::Float64)
    """
    X-polarized.

    Point vector in meters, phase in radians, lambda and length in meters.
    The real part of equation 3 from
    Design Criteria for Near-Field-Focused Planar Arrays,
    [Buffi 2012]
    """
    Rhat = point - element_position

    unit_Rhat = Rhat / norm(Rhat)
    unit_Dxhat = [dipole_length]

    E = (Z_0 / (4 * lambda * R)) * dipole_length * excite_current * sin((2*pi*R/lambda) + phase_shift) * () *

    return E
end

println("result 1: ")

@test dipole_E_field_at_point() == 9

#https://github.com/rdeits/MeshCat.jl#contours
# using MeshCat
# vis = Visualizer()
# open(vis)
# meshcat?
# Visualize a mesh from the level set of a function
# using Meshing
# f = x -> sum(sin, 5 * x)
# sdf = SignedDistanceField(f, HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2)))
# mesh = HomogenousMesh(sdf, MarchingTetrahedra())
# setobject!(vis, mesh,
#            MeshPhongMaterial(color=RGBA{Float32}(1, 0, 0, 0.5)))

# https://github.com/JuliaPlots/Makie.jl
# using LinearAlgebra
#
# using AbstractPlotting
#
# function test(x, y, z)
#     xy = [x, y, z]
#     ((xy') * Matrix(I, 3, 3) * xy) / 20
# end
# x = range(-2pi, stop = 2pi, length = 100)
# scene = Scene()
# # c[4] == fourth argument of the above plotting command
# c = contour!(scene, x, x, x, test, levels = 6, alpha = 0.3, transparency = true)[end]
# xm, ym, zm = minimum(scene_limits(scene))
# contour!(scene, x, x, map(v-> v[1, :, :], c[4]), transformation = (:xy, zm), linewidth = 2)
# heatmap!(scene, x, x, map(v-> v[:, 1, :], c[4]), transformation = (:xz, ym))
# contour!(scene, x, x, map(v-> v[:, :, 1], c[4]), fillrange = true, transformation = (:yz, xm))
# # reorder plots for transparency
# scene.plots[:] = scene.plots[[1, 3, 4, 5, 2]]
# scene

ext: jl


text:
using Printf
using Random;

function sor(n::Int64, A::Array, b::Array, omega::Float64)
	epsilon = 1e-8;
	iter = 0;
	is_convergence = false;
	x = zeros(Float64, n);
	err = zeros(Float64, n);
	while !is_convergence
		is_convergence = true;
		for i = 1:n
			pre_x = x[i];;
			A, b = pivot(n, i, A, b)
			tmp = 0.0;
			s = 0.0;
			for j = 1:n
				if i != j
					s += A[i, j] * x[j]
				end
			end
			tmp = (b[i] - s) / A[i, i];
			x[i] = omega * tmp + (1 - omega) * x[i];
			err[i] = abs(x[i] - pre_x);
		end
		for i = 1:n
			if err[i] > epsilon && iter < 100000
				is_convergence = false
				break;
			end
		end
		iter += 1;
	end
	@printf("iter = %d\n", iter);
	return x;
end

function pivot(size, row, A, b)
    dc = abs(A[row, row]);
    max_row = row
    for i = row:size
        if dc < abs(A[i, row])
            max_row = i;
            dc = abs(A[i, row])
        end
    end
    A[row, :], A[max_row, :] = A[max_row, :], A[row, :];
    b[row], b[max_row] = b[max_row], b[row];
    return A, b
end

function main()
	n = 4;
	omega = 1.0;
	A = rand(Int8, (n, n));
	b = rand(Int16, n);
	for j = 1:n
		for i = 1:n
			if A[j, i] < 0
				@printf("%d x%d ", A[j, i], i);
			else
				@printf("+ %d x%d ", A[j, i], i);
			end
		end
		@printf(" = %d\n", b[j]);
	end
	@time x = sor(n, A, b, omega);
	println();
	for i = 1:n
		@printf("\t\t\tx%d = %.7f\n", i, x[i]);
	end
 end

 main()


text:
mutable struct Compressor{T,P,C,N}
    offsets::Array{Int64,1}
    lengths::Array{Int64,1}
    n::NTuple{N,Int}
    issubarray::Bool
    volume::Array{P,N}
    compressed_volume::Array{UInt32,1}
    compressor::CvxCompressor{N}
end

function Compressor(::Type{T}, ::Type{P}, ::Type{C}, n::NTuple{N,Int}, b::NTuple{N,Int}, scale, nfields, issubarray) where {T,C,P,N}
    volume = zeros(P,ntuple(i->0,N))
    compressed_volume = zeros(UInt32,0)
    compressor = CvxCompressor(b, scale)
    Compressor{T,P,C,N}(zeros(Int64, nfields), zeros(Int64, nfields), n, issubarray, volume, compressed_volume, compressor)
end

function Base.open(compressor::Compressor{T,P,C,N}) where {T,P,C,N}
    compressor.volume = ((T == Float64 && C != Float64) || compressor.issubarray) ? zeros(P,compressor.n) : zeros(P,ntuple(i->0,N))
    compressor.compressed_volume = C == UInt32 ? zeros(UInt32,2*prod(compressor.n)) : zeros(UInt32,0)
    nothing
end

function Base.close(compressor::Compressor{T,P,C,N}) where {T,P,C,N}
    compressor.volume = zeros(T,ntuple(i->0,N))
    compressor.compressed_volume = zeros(UInt32,0)
    nothing
end

copy(compressor::Compressor{T,P,C,N}) where {T,P,C,N} = Compressor{T,P,C,N}(copy(compressor.offsets), copy(compressor.lengths), compressor.n,
    compressor.issubarray, copy(compressor.volume), copy(compressor.compressed_volume), copy(compressor.compressor))

@noinline function compressedwrite_exception(io, nbytes_written, nbytes)
    dir = join(split(split(io.name, " ")[2],"/")[1:end-1], "/")
    dir = dir == "" ? "." : dir
    throw(ErrorException("wrote $(nbytes_written) of $(nbytes) on $(gethostname()) for $(io.name), bytes_available=\n$(run(`df -h $(dir)`))\n"))
end

compressedwrite(io::IO, compressor::Compressor{Float32,Float32,UInt32,N}, ifield::Integer, field::DenseArray{Float32,N}) where {N} = compressedwrite_helper(io, compressor, ifield, field)

function compressedwrite(io::IO, compressor::Compressor{Float32,Float32,UInt32,N}, ifield::Integer, field::DenseArray{Float32,N}, rng::NTuple{N,UnitRange{Int}}) where N
    copyto!(compressor.volume, view(field, rng...))
    compressedwrite_helper(io, compressor, ifield, compressor.volume)
end

function compressedwrite(io::IO, compressor::Compressor{Float64,Float32,UInt32,N}, ifield::Integer, field::DenseArray{Float64,N}) where N
    copyto!(compressor.volume, field)
    compressedwrite_helper(io, compressor, ifield, compressor.volume)
end

function compressedwrite(io::IO, compressor::Compressor{Float64,Float32,UInt32,N}, ifield::Integer, field::DenseArray{Float64,N}, rng::NTuple{N,UnitRange{Int}}) where N
    copyto!(compressor.volume, view(field, rng...))
    compressedwrite_helper(io, compressor, ifield, compressor.volume)
end

function compressedwrite_helper(io::IO, compressor::Compressor{T,P,UInt32,N}, ifield::Integer, field::DenseArray{Float32,N}) where {T,P,N}
    compressor.lengths[ifield] = compress!(compressor.compressed_volume, compressor.compressor, field)
    nbytes = unsafe_write(io, convert(Ptr{UInt8},pointer(compressor.compressed_volume)), compressor.lengths[ifield])
    nbytes == compressor.lengths[ifield] || compressedwrite_exception(io, nbytes, compressor.lengths[ifield])
    if ifield < length(compressor.offsets)
        @inbounds compressor.offsets[ifield+1] = compressor.offsets[ifield] + compressor.lengths[ifield]
    end
end

function compressedwrite(io::IO, compressor::Compressor{T,T,T,N}, ifield::Integer, field::DenseArray{T,N}) where {T<:Union{Float32,Float64},N}
    nbytes = write(io, field)
    nbytes == length(field)*sizeof(T) || compressedwrite_exception(io, nbytes, compressor.lengths[ifield])
end

function compressedwrite(io::IO, compressor::Compressor{T,T,T,N}, ifield::Integer, field::DenseArray{T,N}, rng::NTuple{N,UnitRange{Int}}) where {T<:Union{Float32,Float64},N}
    _field = view(field, rng...)
    copyto!(compressor.volume, _field)
    nbytes = write(io, compressor.volume)
    nbytes == length(_field)*sizeof(T) || compressedwrite_exception(io, nbytes, compressor.lengths[ifield])
end

function compressedwrite(io::IO, compressor::Compressor{Float64,Float32,Float32,N}, ifield::Integer, field::DenseArray{Float64,N}) where N
    copyto!(compressor.volume, field)
    nbytes = write(io, compressor.volume)
    nbytes == length(field)*sizeof(Float32) || compressedwrite_exception(io, nbytes, compressor.lengths[ifield])
end

function compressedwrite(io::IO, compressor::Compressor{Float64,Float32,Float32,N}, ifield::Integer, field::DenseArray{Float64,N}, rng::NTuple{N,UnitRange{Int}}) where N
    _field = view(field, rng...)
    copyto!(compressor.volume, _field)
    nbytes = write(io, compressor.volume)
    nbytes == length(_field)*sizeof(Float32) || compressedwrite_exception(io, nbytes, compressor.lengths[ifield])
end

compressedread!(io::IO, compressor::Compressor{Float32,Float32,UInt32,N}, ifield::Integer, field::DenseArray{Float32,N}) where {N} = compressedread_helper!(io, compressor, ifield, field)

function compressedread!(io::IO, compressor::Compressor{Float32,Float32,UInt32,N}, ifield::Integer, field::DenseArray{Float32,N}, rng::NTuple{N,UnitRange{Int}}) where N
    compressedread_helper!(io, compressor, ifield, compressor.volume)
    copyto!(view(field,rng...), compressor.volume)
end

function compressedread!(io::IO, compressor::Compressor{Float64,Float32,UInt32,N}, ifield::Integer, field::DenseArray{Float64,N}) where N
    compressedread_helper!(io, compressor, ifield, compressor.volume)
    copyto!(field, compressor.volume)
end

function compressedread!(io::IO, compressor::Compressor{Float64,Float32,UInt32,N}, ifield::Integer, field::DenseArray{Float64,N}, rng::NTuple{N,UnitRange{Int}}) where N
    compressedread_helper!(io, compressor, ifield, compressor.volume)
    copyto!(view(field, rng...), compressor.volume)
end

function compressedread_helper!(io::IO, compressor::Compressor{T,Float32,UInt32,N}, ifield::Integer, field::DenseArray{Float32,N}) where {T,N}
    seek(io, compressor.offsets[ifield])
    compressed_volume_ptr = convert(Ptr{UInt8}, pointer(compressor.compressed_volume))
    unsafe_read(io, compressed_volume_ptr, compressor.lengths[ifield])
    decompress!(field, compressor.compressor, compressor.compressed_volume, compressor.lengths[ifield])
end

function compressedread!(io::IO, compressor::Compressor{T,T,T,N}, ifield::Integer, field::DenseArray{T,N}) where {T<:Union{Float32,Float64},N}
    seek(io, (ifield-1)*length(field)*sizeof(T))
    read!(io, field)
end

function compressedread!(io::IO, compressor::Compressor{T,T,T,N}, ifield::Integer, field::DenseArray{T,N}, rng::NTuple{N,UnitRange{Int}}) where {T<:Union{Float32,Float64},N}
    _field = view(field, rng...)
    seek(io, (ifield-1)*length(_field)*sizeof(T))
    read!(io, compressor.volume)
    copyto!(_field, compressor.volume)
end

function compressedread!(io::IO, compressor::Compressor{Float64,Float32,Float32,N}, ifield::Integer, field::DenseArray{Float64,N}) where N
    seek(io, (ifield-1)*length(field)*sizeof(Float32))
    read!(io, compressor.volume)
    copyto!(field, compressor.volume)
end

function compressedread!(io::IO, compressor::Compressor{Float64,Float32,Float32,N}, ifield::Integer, field::DenseArray{Float64,N}, rng::NTuple{N,UnitRange{Int}}) where N
    _field = view(field, rng...)
    seek(io, (ifield-1)*length(_field)*sizeof(Float32))
    read!(io, compressor.volume)
    copyto!(_field, compressor.volume)
end

function comptype(comptype, T)
    comptype = comptype == nothing ? T : comptype
    if T == Float64
        @assert comptype == Float64 || comptype == Float32 || comptype == UInt32
    else
        @assert comptype == Float32 || comptype == UInt32
    end
    P = (T == Float64 && comptype == Float64) ? Float64 : Float32
    comptype, P
end

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
writeline(file, s) = write(file, string(s, "\r\n"))
for i = 1:30
	filename = "ambulances_$i.csv"
	open(filename, "w") do file
		writeline(file, "ambulances,,")
		writeline(file, "index,stationIndex,class")
		for j = 1:i
			writeline(file, "$j,1,1")
		end
	end
end


text:
""" get_obs_loc(obs::Symbol)
 Returns a Dict with long & lat (degrees) and elevation (km)
 Warning: Currently only has info for :WIYN and :HARPSN.
"""
function get_obs_loc(obs::Symbol)
	valid_obs = [:WIYN, :HARPSN]
	@assert obs ∈ valid_obs
	if obs == :HARPSN
		loc = Dict("lon"=> -17.88905, "lat"=> 28.754, "elevation"=> 2.3872)
	elseif obs == :WIYN
		loc = Dict("lon"=> -111.600562, "lat"=> 31.958092, "elevation"=> 2.091)
	else
		@error("Don't have coordinates for obs = " * string(obs))
	end
	return loc
end

function calc_solar_alt(time::Real; obs::Symbol)
	loc = get_obs_loc(obs)
  #=
	obs = AstroLib.observatories["kpno"]
	lat = obs.latitude
	long = obs.longitude
	=#
  #=
	 Values for specifically for WIYN
  long = -111.600562
  lat =  31.958092
  alt =  2091.0
	=#
  (ra_sun, dec_sun, ) = AstroLib.sunpos(time)
  lst = AstroLib.ct2lst(loc["lon"],time)
  (alt, az) = AstroLib.hadec2altaz(lst-ra_sun/360*24, dec_sun,loc["lat"])
  return alt
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct SimpleFullMachine <: Machine
        R::Float64
        R_f::Float64
        R_1d::Float64
        R_1q::Float64
        L_d::Float64
        L_q::Float64
        L_ad::Float64
        L_aq::Float64
        L_f1d::Float64
        L_ff::Float64
        L_1d::Float64
        L_1q::Float64
        MVABase::Float64
        ext::Dict{String, Any}
        inv_d_fluxlink::Array{Float64,2}
        inv_q_fluxlink::Array{Float64,2}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameter of a full order flux stator-rotor model without zero sequence flux in the stator.
 The derivative of stator fluxes (ψd and ψq) is neglected. This is standard when
 transmission network dynamics is neglected. Only one q-axis damping circuit
 is considered. All per unit are in machine per unit.
 Refer to Chapter 3 of Power System Stability and Control by P. Kundur or Chapter 11 of Power System Dynamics: Stability and Control, by J. Machowski, J. Bialek and J. Bumby, for more details.
 Note that the models are somewhat different (but equivalent) due to the different Park Transformation used in both books.

# Arguments
- `R::Float64`: Resistance after EMF in machine per unit, validation range: (0, nothing)
- `R_f::Float64`: Field rotor winding resistance in per unit, validation range: (0, nothing)
- `R_1d::Float64`:  Damping rotor winding resistance on d-axis in per unit. This value is denoted as RD in Machowski., validation range: (0, nothing)
- `R_1q::Float64`: Damping rotor winding resistance on q-axis in per unit. This value is denoted as RQ in Machowski., validation range: (0, nothing)
- `L_d::Float64`: Inductance of fictitious damping that represent the effect of the three-phase stator winding in the d-axis of the rotor, in per unit. This value is denoted as L_ad + L_l in Kundur (and Ld in Machowski)., validation range: (0, nothing)
- `L_q::Float64`: Inductance of fictitious damping that represent the effect of the three-phase stator winding in the q-axis of the rotor, in per unit. This value is denoted as L_aq + L_l in Kundur., validation range: (0, nothing)
- `L_ad::Float64`: Mutual inductance between stator winding and rotor field (and damping) winding inductance on d-axis, in per unit, validation range: (0, nothing)
- `L_aq::Float64`: Mutual inductance between stator winding and rotor damping winding inductance on q-axis, in per unit, validation range: (0, nothing)
- `L_f1d::Float64`: Mutual inductance between rotor field winding and rotor damping winding inductance on d-axis, in per unit, validation range: (0, nothing)
- `L_ff::Float64`: Field rotor winding inductance, in per unit, validation range: (0, nothing)
- `L_1d::Float64`: Inductance of the d-axis rotor damping circuit, in per unit, validation range: (0, nothing)
- `L_1q::Float64`: Inductance of the q-axis rotor damping circuit, in per unit, validation range: (0, nothing)
- `MVABase::Float64`: Nominal Capacity in MVA, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `inv_d_fluxlink::Array{Float64,2}`: Equations 3.127, 3.130, 3.131 From Kundur
- `inv_q_fluxlink::Array{Float64,2}`: Equations 3.128, 3.132 From Kundur
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct SimpleFullMachine <: Machine
    "Resistance after EMF in machine per unit"
    R::Float64
    "Field rotor winding resistance in per unit"
    R_f::Float64
    " Damping rotor winding resistance on d-axis in per unit. This value is denoted as RD in Machowski."
    R_1d::Float64
    "Damping rotor winding resistance on q-axis in per unit. This value is denoted as RQ in Machowski."
    R_1q::Float64
    "Inductance of fictitious damping that represent the effect of the three-phase stator winding in the d-axis of the rotor, in per unit. This value is denoted as L_ad + L_l in Kundur (and Ld in Machowski)."
    L_d::Float64
    "Inductance of fictitious damping that represent the effect of the three-phase stator winding in the q-axis of the rotor, in per unit. This value is denoted as L_aq + L_l in Kundur."
    L_q::Float64
    "Mutual inductance between stator winding and rotor field (and damping) winding inductance on d-axis, in per unit"
    L_ad::Float64
    "Mutual inductance between stator winding and rotor damping winding inductance on q-axis, in per unit"
    L_aq::Float64
    "Mutual inductance between rotor field winding and rotor damping winding inductance on d-axis, in per unit"
    L_f1d::Float64
    "Field rotor winding inductance, in per unit"
    L_ff::Float64
    "Inductance of the d-axis rotor damping circuit, in per unit"
    L_1d::Float64
    "Inductance of the q-axis rotor damping circuit, in per unit"
    L_1q::Float64
    "Nominal Capacity in MVA"
    MVABase::Float64
    ext::Dict{String, Any}
    "Equations 3.127, 3.130, 3.131 From Kundur"
    inv_d_fluxlink::Array{Float64,2}
    "Equations 3.128, 3.132 From Kundur"
    inv_q_fluxlink::Array{Float64,2}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function SimpleFullMachine(R, R_f, R_1d, R_1q, L_d, L_q, L_ad, L_aq, L_f1d, L_ff, L_1d, L_1q, MVABase, ext=Dict{String, Any}(), )
    SimpleFullMachine(R, R_f, R_1d, R_1q, L_d, L_q, L_ad, L_aq, L_f1d, L_ff, L_1d, L_1q, MVABase, ext, inv([[-L_d L_ad L_ad]; [-L_ad L_ff L_f1d]; [-L_ad L_f1d L_1d]]), inv([[-L_q L_aq]; [-L_aq L_1q]]), [:ψf, :ψ1d, :ψ1q], 3, InfrastructureSystemsInternal(), )
end

function SimpleFullMachine(; R, R_f, R_1d, R_1q, L_d, L_q, L_ad, L_aq, L_f1d, L_ff, L_1d, L_1q, MVABase, ext=Dict{String, Any}(), )
    SimpleFullMachine(R, R_f, R_1d, R_1q, L_d, L_q, L_ad, L_aq, L_f1d, L_ff, L_1d, L_1q, MVABase, ext, )
end

# Constructor for demo purposes; non-functional.
function SimpleFullMachine(::Nothing)
    SimpleFullMachine(;
        R=0,
        R_f=0,
        R_1d=0,
        R_1q=0,
        L_d=1.0,
        L_q=1.0,
        L_ad=2.0,
        L_aq=1.0,
        L_f1d=1.0,
        L_ff=2.0,
        L_1d=1.0,
        L_1q=2.0,
        MVABase=0,
        ext=Dict{String, Any}(),
    )
end

"""Get SimpleFullMachine R."""
get_R(value::SimpleFullMachine) = value.R
"""Get SimpleFullMachine R_f."""
get_R_f(value::SimpleFullMachine) = value.R_f
"""Get SimpleFullMachine R_1d."""
get_R_1d(value::SimpleFullMachine) = value.R_1d
"""Get SimpleFullMachine R_1q."""
get_R_1q(value::SimpleFullMachine) = value.R_1q
"""Get SimpleFullMachine L_d."""
get_L_d(value::SimpleFullMachine) = value.L_d
"""Get SimpleFullMachine L_q."""
get_L_q(value::SimpleFullMachine) = value.L_q
"""Get SimpleFullMachine L_ad."""
get_L_ad(value::SimpleFullMachine) = value.L_ad
"""Get SimpleFullMachine L_aq."""
get_L_aq(value::SimpleFullMachine) = value.L_aq
"""Get SimpleFullMachine L_f1d."""
get_L_f1d(value::SimpleFullMachine) = value.L_f1d
"""Get SimpleFullMachine L_ff."""
get_L_ff(value::SimpleFullMachine) = value.L_ff
"""Get SimpleFullMachine L_1d."""
get_L_1d(value::SimpleFullMachine) = value.L_1d
"""Get SimpleFullMachine L_1q."""
get_L_1q(value::SimpleFullMachine) = value.L_1q
"""Get SimpleFullMachine MVABase."""
get_MVABase(value::SimpleFullMachine) = value.MVABase
"""Get SimpleFullMachine ext."""
get_ext(value::SimpleFullMachine) = value.ext
"""Get SimpleFullMachine inv_d_fluxlink."""
get_inv_d_fluxlink(value::SimpleFullMachine) = value.inv_d_fluxlink
"""Get SimpleFullMachine inv_q_fluxlink."""
get_inv_q_fluxlink(value::SimpleFullMachine) = value.inv_q_fluxlink
"""Get SimpleFullMachine states."""
get_states(value::SimpleFullMachine) = value.states
"""Get SimpleFullMachine n_states."""
get_n_states(value::SimpleFullMachine) = value.n_states
"""Get SimpleFullMachine internal."""
get_internal(value::SimpleFullMachine) = value.internal


text:
function airyai(x::ArbReal{P}, prec::Int=P) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_hypgeom_airy), Cvoid,
              (Ref{ArbReal}, Ref{Cvoid}, Ref{Cvoid}, Ref{Cvoid}, Ref{ArbReal}, Clong),
              z, C_NULL, C_NULL, C_NULL, x, prec)
    return z
end

function airyaiprime(x::ArbReal{P}, prec::Int=P) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_hypgeom_airy), Cvoid,
              (Ref{Cvoid}, Ref{ArbReal}, Ref{Cvoid}, Ref{Cvoid}, Ref{ArbReal}, Clong),
              C_NULL, z, C_NULL, C_NULL, x, prec)
    return z
end

function airybi(x::ArbReal{P}, prec::Int=P) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_hypgeom_airy), Cvoid,
              (Ref{Cvoid}, Ref{Cvoid}, Ref{ArbReal}, Ref{Cvoid}, Ref{ArbReal}, Clong),
              C_NULL, C_NULL, z, C_NULL, x, prec)
    return z
end

function airybiprime(x::ArbReal{P}, prec::Int=P) where {P}
    z = ArbReal{P}()
    ccall(@libarb(arb_hypgeom_airy), Cvoid,
              (Ref{Cvoid}, Ref{Cvoid}, Ref{Cvoid}, Ref{ArbReal}, Ref{ArbReal}, Clong),
              C_NULL, C_NULL, C_NULL, z, x, prec)
    return z
end

for F in (:airyai, :airyaiprime, :airybi, :airybiprime)
    @eval begin
        function $F(x::ArbFloat{P}, prec::Int=P) where {P}
            r = ArbReal{P}(x)
            r = $F(r, prec)
            z = ArbFloat{P}(r)
            return r
        end
    end
end

function airyai(x::ArbComplex{P}, prec::Int=P) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_hypgeom_airy), Cvoid,
              (Ref{ArbComplex}, Ref{Cvoid}, Ref{Cvoid}, Ref{Cvoid}, Ref{ArbComplex}, Clong),
              z, C_NULL, C_NULL, C_NULL, x, prec)
    return z
end

function airyaiprime(x::ArbComplex{P}, prec::Int=P) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_hypgeom_airy), Cvoid,
              (Ref{Cvoid}, Ref{ArbComplex}, Ref{Cvoid}, Ref{Cvoid}, Ref{ArbComplex}, Clong),
              C_NULL, z, C_NULL, C_NULL, x, prec)
    return z
end

function airybi(x::ArbComplex{P}, prec::Int=P) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_hypgeom_airy), Cvoid,
              (Ref{Cvoid}, Ref{Cvoid}, Ref{ArbComplex}, Ref{Cvoid}, Ref{ArbComplex}, Clong),
              C_NULL, C_NULL, z, C_NULL, x, prec)
    return z
end

function airybiprime(x::ArbComplex{P}, prec::Int=P) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_hypgeom_airy), Cvoid,
              (Ref{Cvoid}, Ref{Cvoid}, Ref{Cvoid}, Ref{ArbComplex}, Ref{ArbComplex}, Clong),
              C_NULL, C_NULL, C_NULL, z, x, prec)
    return z
end


text:
using ITensors
using ITensorVisualization
using LayeredLayouts
using Graphs

N = 10
layers = 10
ndelete = 0

s = siteinds("Qubit", N)
layer(N, start) = [("CX", i, i + 1) for i in start:2:(N - 1)]
layer(N) = append!(layer(N, 1), layer(N, 2))
layer_N = layer(N)
gates = []
for _ in 1:layers
  append!(gates, layer_N)
end

for _ in 1:ndelete
  deleteat!(gates, rand(eachindex(gates)))
end

U, s̃ = circuit_network(gates, s)
ψ = prod(MPS(s))
ψ̃ = prod(MPS(s̃))
tn = [ψ, U..., ψ̃]

edge_labels = (; plevs=true)
layout(g) = layered_layout(solve_positions(Zarate(), g))
@visualize fig tn arrow_show = true edge_labels = edge_labels layout = layout

fig


text:
import CImGui.LibCImGui:
    GetPlotMousePos,
    GetPlotLimits,
    GetPlotQuery,
    GetLastItemColor,
    NextColormapColor

function SetNextPlotTicksX(values::Vector{<:Real}, n_ticks::Integer,
                           labels = [""], show_default::Bool = false)

    eltype(values) !== Float64 && (values = Float64.(values))

    ImPlot_SetNextPlotTicksXdoublePtr(values, Cint(n_ticks), labels, show_default)
end

function SetNextPlotTicksX(x_min, x_max, n_ticks::Integer,
                           labels = [""], show_default::Bool = false)

    typeof(x_min) !== Float64 && (x_min = Float64(x_min))
    typeof(x_max) !== Float64 && (x_max = Float64(x_max))

    ImPlot_SetNextPlotTicksXdouble(x_min, x_max, Cint(n_ticks), labels, show_default)
end

function SetNextPlotTicksY(values::Vector{<:Real}, n_ticks::Integer,
                           labels = [""], show_default::Bool = false,
                           y_axis::Integer = 0)

    eltype(values) !== Float64 && (values = Float64.(values))

    ImPlot_SetNextPlotTicksYdoublePtr(values, Cint(n_ticks), labels, show_default,
                                          Cint(y_axis))
end

function SetNextPlotTicksY(y_min, y_max, n_ticks::Integer,
                           labels = [""], show_default::Bool = false,
                           y_axis::Integer = 0)

    typeof(y_min) !== Float64 && (y_min = Float64(y_min))
    typeof(y_max) !== Float64 && (y_max = Float64(y_max))

    ImPlot_SetNextPlotTicksYdouble(y_min, y_max, Cint(n_ticks), labels, show_default,
                                       Cint(y_axis))
end

function GetPlotPos()
    out = Ref{ImVec2}()
    ImPlot_GetPlotPos(out)
    return out[]
end

function GetPlotSize()
    out = Ref{ImVec2}()
    ImPlot_GetPlotSize(out)
    return out[]
end

SetLegendLocation(location, orientation) = ImPlot_SetLegendLocation(location, orientation, false)
SetLegendLocation(location) = ImPlot_SetLegendLocation(location, ImPlot_ImPlotOrientation_Vertical, false)

# ImPlot_SetNextPlotLimits(xmin, xmax, ymin, ymax) = ImPlot_SetNextPlotLimits(xmin, xmax, ymin, ymax, ImGuiCond_Once)
SetNextPlotLimits(xmin, xmax, ymin, ymax, cond=ImGuiCond_Once) = ImPlot_SetNextPlotLimits(xmin, xmax, ymin, ymax, cond)

function PixelsToPlot(pix::ImVec2, y_axis = IMPLOT_AUTO)
    out = Ref(ImPlotPoint())
    ImPlot_PixelsToPlotVec2(out, pix, y_axis)
    return out[]
end

function PixelsToPlot(x::Real, y::Real, y_axis = IMPLOT_AUTO)
    out = Ref(ImPlotPoint())
    ImPlot_PixelsToPlotFloat(out, x, y, y_axis)
    return out[]
end
function PlotToPixels(plt::ImPlotPoint, y_axis = IMPLOT_AUTO)
    out = Ref{ImVec2}()
    ImPlot_PlotToPixelsPlotPoInt(out, plt, y_axis)
    return out[]
end
function PlotToPixels(x::Real, y::Real, y_axis = IMPLOT_AUTO)
    out = Ref{ImVec2}()
    ImPlot_PlotToPixelsdouble(out, x, y, y_axis)
    return out[]
end

function ImPlot_GetPlotMousePos(y_axis = IMPLOT_AUTO)
    out = Ref(ImPlotPoint())
    ImPlot_GetPlotMousePos(out, y_axis)
    return out[]
end

function ImPlot_GetPlotLimits(y_axis = IMPLOT_AUTO)
    out = Ref(ImPlotLimits())
    ImPlot_GetPlotLimits(out, y_axis)
    return out[]
end

function ImPlot_GetPlotQuery(y_axis = IMPLOT_AUTO)
    out = Ref(ImPlotLimits())
    ImPlot_GetPlotQuery(out, y_axis)
    return out[]
end

function ImPlot_GetLastItemColor()
    out = Ref{ImVec4}()
    ImPlot_GetLastItemColor(out)
    return out[]
end
function ImPlot_GetColormapColor(index)
    out = Ref{ImVec4}()
    ImPlot_GetColormapColor(out, index)
    return out[]
end
function ImPlot_LerpColormap(t)
    out = Ref{ImVec4}()
    ImPlot_LerpColormap(out, t)
    return out[]
end
function ImPlot_NextColormapColor()
    out = Ref{ImVec4}()
    ImPlot_NextColormapColor(out)
    return out[]
end

# function Contains(range::ImPlot_ImPlotRange, value)
#     return value >= range.Min && value <= range.Max
# end
# function Contains(limits::ImPlot_ImPlotLimits, x, y)
#     return Contains(limits.X, x) && Contains(limits.Y, y)
# end
# function Contains(limits::ImPlot_ImPlotLimits, p::ImPlot_ImPlotPoint)
#     return Contains(limits.X, p.x) && Contains(limits.Y, p.y)
# end

function DragLineX(id::String, x_value, show_label::Bool = true, col::ImVec4 = IMPLOT_AUTO_COL, thickness::Real = 1)
    ImPlot_DragLineX(id, x_value, show_label, col, thickness)
end

function DragLineY(id::String, y_value, show_label::Bool = true, col::ImVec4 = IMPLOT_AUTO_COL, thickness::Real = 1)
    ImPlot_DragLineY(id, y_value, show_label, col, thickness)
end

function DragPoint(id::String, x, y, show_label::Bool = true, col::ImVec4 = IMPLOT_AUTO_COL, radius::Real = 4)
    ImPlot_DragPoint(id, x, y, show_label, col, radius)
end

SetNextLineStyle(col = IMPLOT_AUTO_COL, weight = IMPLOT_AUTO)  = ImPlot_SetNextLineStyle(col,weight)
SetNextFillStyle(col = IMPLOT_AUTO_COL, alpha_mod = IMPLOT_AUTO) = ImPlot_SetNextFillStyle(col, alpha_mod)
SetNextMarkerStyle(marker = IMPLOT_AUTO, size = IMPLOT_AUTO, fill = IMPLOT_AUTO_COL, weight = IMPLOT_AUTO, outline = IMPLOT_AUTO_COL) = ImPlot_SetNextMarkerStyle(marker, size, fill, weight, outline)
SetNextErrorBarStyle(col = IMPLOT_AUTO_COL, size = IMPLOT_AUTO, weight = IMPLOT_AUTO) = ImPlot_SetNextErrorBarStyle(col, size, weight)


text:
# Tree API implementations for builtin types

# Types which are iterable but shouldn't be considered tree-iterable
children(x::Number) = ()
children(x::Char) = ()
children(x::Task) = ()
children(x::AbstractString) = ()

# Expr
children(x::Expr) = x.args

function printnode(io::IO, x::Expr)
    print(io, "Expr(")
    show(io, x.head)
    print(io, ")")
end

# AbstractDict
printnode(io::IO, kv::Pair{K,V}) where {K,V} = printnode(io,kv[1])
children(kv::Pair{K,V}) where {K,V} = (kv[2],)

# For potentially-large containers, just show the type
printnode(io::IO, ::T) where T <: Union{AbstractArray, AbstractDict} = print(io, T)

if isdefined(Core.Compiler, :Timings)
    children(t::Core.Compiler.Timings.Timing) = t.children
    printnode(io::IO, t::Core.Compiler.Timings.Timing) = print(io, t.time/10^6, "ms: ", t.mi_info)
    nodetype(t::Core.Compiler.Timings.Timing) = Core.Compiler.Timings.Timing
end


text:
using Statistics

X = rand(50, 20)

# True generating function has intercept, 2 factors and a lag range of 3-5.
genfn(r, X) = 2.56 - 7.8*mean(X[(r-5):(r-4), 2]) + 9.3*X[r-3, 11] + (0.01 * rand())

y = Float64[genfn(r, X) for r in 6:size(X, 1)]

# Now in reality we don't know the real model se we will have to train it
# from data. Let's do it and try to recover an interpretable model.
using InterpretableModels

# We don't know how many factors so let's aim for 3 and a lag range of 2-10
m = MovingAverageLITRM(3, size(X, 2); minlag=2, maxlag=10)

fit!(m, X, y; MaxTime = 10.0)

println(m) # The model should have been recovered quite well...


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.1
#     language: julia
#     name: julia-1.6
# ---

# %%
InteractiveUtils.code_llvm(@nospecialize(f), @nospecialize(types=Tuple); raw=false, dump_module=false, optimize=true, debuginfo::Symbol=:default) = sprint() do io
    code_llvm(io, f, types; raw=raw, dump_module=dump_module, optimize=optimize, debuginfo=debuginfo)
end

# %%
s = @code_llvm debuginfo=:none Base.Math.sin_kernel(0.5)

# %%
s |> typeof

# %%
print(s)

# %%

max_stars_repo_path: 0011/Overwrite code_llvm to return a string.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using ITensors
using ITensorTDVP

include("05_utils.jl")

function heisenberg(N)
  os = OpSum()
  for j in 1:(N - 1)
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
    os += "Sz", j, "Sz", j + 1
  end
  return os
end

N = 10
cutoff = 1e-12
outputlevel = 1
nsteps = 10
time_steps = [n ≤ 2 ? -0.2im : -0.1im for n in 1:nsteps]

obs = Observer("times" => (; current_time) -> current_time, "psis" => (; psi) -> psi)

s = siteinds("S=1/2", N; conserve_qns=true)
H = MPO(heisenberg(N), s)

psi0 = MPS(s, n -> isodd(n) ? "Up" : "Dn")
psi = tdvp_nonuniform_timesteps(
  ProjMPO(H), psi0; time_steps, cutoff, outputlevel, (step_observer!)=obs
)

res = results(obs)
times = res["times"]
psis = res["psis"]

println("\nResults")
println("=======")
print("step = ", 0)
print(", time = ", zero(ComplexF64))
print(", ⟨Sᶻ⟩ = ", round(expect(psi0, "Sz"; sites=N ÷ 2); digits=3))
println()
for n in 1:length(times)
  print("step = ", n)
  print(", time = ", round(times[n]; digits=3))
  print(", ⟨Sᶻ⟩ = ", round(expect(psis[n], "Sz"; sites=N ÷ 2); digits=3))
  println()
end


text:
"""
ptcsol(F!, x0, FS, FPS, J! = diffjac!; rtol=1.e-6, atol=1.e-12,
             maxit=20, delta0=1.e-6, dx=1.e-7, pdata = nothing, jfact = klfact,
               printerr = true, keepsolhist = false, jknowsdt = false)

C. T. Kelley, 2021

Julia versions of the nonlinear solvers from my SIAM books.
Herewith: some new stuff ==> ptcsol

PTC finds the steady-state solution of u' = -F(u), u(0) = u_0.
The - sign is a convention.

You must allocate storage for the function and Jacobian in advance
--> in the calling program <-- ie. in FS and FPS

Inputs:\n
- F!: function evaluation, the ! indicates that F! overwrites FS, your
    preallocated storage for the function.\n
    So, FV=F!(FV,x) or FV=F!(FV,x,pdata) returns FV=F(x)

- x0: initial iterate\n

- FS: Preallocated storage for function. It is an N x 1 column vector.\n
You may dimension it as (n,) or (n,1). (n,) is best, but the
solvers can deal with it either way.

- FPS: preallocated storage for Jacobian. It is an N x N matrix\n
  If FPS is sparse, you __must__ allocate storage for the diagonal so
  I will have room to put 1/dt in there.

- J!: Jacobian evaluation, the ! indicates that J! overwrites FPS, your
    preallocated storage for the Jacobian. If you leave this out the
    default is a finite difference Jacobian.\n
    So, FP=J!(FP,FV,x) or FP=J!(FP,FV,x,pdata) returns FP=F'(x);
    (FP,FV, x) must be the argument list, even if FP does not need FV.
    One reason for this is that the finite-difference Jacobian
    does and that is the default in the solver.

    You may have a better way to add (1/dt) I to your Jacobian. If you
    want to do this yourself then your Jacobian function should be
    FP=J!(FP,FV,x,dt) or FP=J!(FP,FV,x,dt,pdata) and return
    F'(x) + (1.0/dt)*I. \n
    You will also have to set the kwarg __jknowsdt__ to true.

- Precision: Lemme tell ya 'bout precision. I designed this code for
    full precision
    functions and linear algebra in any precision you want. You can declare
    FPS as Float64, Float32, or Float16 and ptcsol will do the right thing if
    YOU do not destroy the declaration in your J! function. I'm amazed
    that this works so easily. If the Jacobian is reasonably well
    conditioned, you can cut the cost of Jacobian factorization and
    storage in half with no loss. For large dense Jacobians and inexpensive
    functions, this is a good deal.\n
    BUT ... There is very limited support for direct sparse solvers in
    anything other than Float64. I recommend that you only use Float64
    with direct sparse solvers unless you really know what you're doing. I
    have a couple examples in the notebook, but watch out.

----------------------

Keyword Arguments (kwargs):\n

rtol and atol: relative and absolute error tolerances\n

delta0: initial pseudo time step. The default value of 1.e-3 is a bit conservative
and is one option you really should play with. Look at the example
where I set it to 1.0!\n

maxit: limit on nonlinear iterations, default=100. \n
This is coupled to delta0. If your choice of delta0 is too small (conservative)
then you'll need many iterations to converge and will need a larger
value of maxit

For PTC you'll need more iterations than for a straight-up
nonlinear solve. This is part of the price for finding the
stable solution.

dx: default = 1.e-7\n
difference increment in finite-difference derivatives
      h=dx*norm(x)+1.e-6

pdata:\n
precomputed data for the function/Jacobian.
Things will go better if you use this rather than hide the data
in global variables within the module for your function/Jacobian

jfact: default = klfact (tries to figure out best choice) \n
If your Jacobian has any special structure, please set jfact
to the correct choice for a factorization.

I use jfact when I call PTCUpdate to evaluate the
Jacobian (using your J!) and factor it. The default is to use
klfact (an internal function) to do something reasonable.
For general matrices, klfact picks lu! to compute an LU factorization
and share storage with the Jacobian.  You may change LU to something else by,
for example, setting jfact = cholseky! if your Jacobian is spd.

klfact knows about banded matrices and picks qr. You should,
however RTFM, allocate the extra two upper bands, and use jfact=qr!
to override klfact.

If you give me something that klfact does not know how to dispatch on,
then nothing happens. I just return the original Jacobian matrix and
ptcsol will use backslash to compute the Newton step.

I know that this is probably not optimal in your situation, so it is
good to pick something else, like jfact = lu.

printerr: default = true\n
I print a helpful message when the solver fails. To suppress that
message set printerr to false.

keepsolhist: default = false\n
Set this to true to get the history of the iteration in the output
tuple. This is on by default for scalar equations and off for systems.
Only turn it on if you have use for the data, which can get REALLY LARGE.

jknowsdt: default = false\n
Set this to true if your Jacobian evaluation function retursn
F'(x) + (1/dt) I. You'll also need to follow the rules above for
the Jacobian evaluation function. I do not recommend this and if
your Jacobian is anything other than a matrix I can't promise
anything. I've tested this for matrix outputs only.

Output:\n
A named tuple (solution, functionval, history, stats, idid,
               errcode, solhist)
where

solution = converged result
functionval = F(solution)
history = the vector of residual norms (||F(x)||) for the iteration

Unlike nsol, nsoli, or even ptcsoli, ptcsol has a fixed cost per
iteration of one function, one Jacobian, and one Factorization. Hence
iteration statistics are not interesting and not in the output.

idid=true if the iteration succeeded and false if not.

errcode = 0 if if the iteration succeeded
        = -1 if the initial iterate satisfies the termination criteria
        = 10 if no convergence after maxit iterations

solhist:\n
This is the entire history of the iteration if you've set
keepsolhist=true\n
solhist is an N x K array where N is the length of x and K is the number
of iteration + 1. So, for scalar equations, it's a row vector.

### Example for ptcsol

#### The buckling beam problem.
You'll need to use TestProblems for this to work.

```jldoctest
julia> using SIAMFANLEquations.TestProblems

julia> n=63; maxit=1000; delta = 0.01; lambda = 20.0;

julia> bdata = beaminit(n, 0.0, lambda);

julia> x = bdata.x; u0 = x .* (1.0 .- x) .* (2.0 .- x); u0 .*= exp.(-10.0 * u0);

julia> FS = copy(u0); FPS = copy(bdata.D2);

julia> pout = ptcsol( FBeam!, u0, FS, FPS, BeamJ!; rtol = 1.e-10, pdata = bdata,
                delta0 = delta, maxit = maxit);

julia> # It takes a few iterations to get there.
       length(pout.history)
25

julia> [pout.history[1:5] pout.history[21:25]]
5×2 Array{Float64,2}:
 6.31230e+01  9.75412e-01
 7.52624e+00  8.35295e-02
 8.31545e+00  6.58797e-04
 3.15455e+01  4.12697e-08
 3.66566e+01  6.75094e-12

julia> # We get the nonnegative stedy state.
       norm(pout.solution,Inf)
2.19086e+00
```

"""
function ptcsol(
    F!,
    x0,
    FS = [],
    FPS = [],
    J! = diffjac!;
    rtol = 1.e-6,
    atol = 1.e-12,
    maxit = 20,
    delta0 = 1.e-6,
    dx = 1.e-7,
    pdata = nothing,
    jfact = klfact,
    printerr = true,
    keepsolhist = false,
    jknowsdt = false,
)
    itc = 0
    idid = true
    #
    #   Initialize the iteration
    #   As with the other codes, ItRules packages all the details of
    #   the problem so it's easy to pass them around.
    #
    (ItRules, x, n, solhist) =
        PTCinit(x0, dx, F!, J!, delta0, maxit, pdata, jfact, keepsolhist, jknowsdt)
    #
    # First Evaluation of the function. Initialize the iteration history.
    # Fix the tolerances for convergence and define the derivative FPF
    # outside of the main loop for scoping.
    #
    FS = EvalF!(F!, FS, x, pdata)
    resnorm = norm(FS)
    tol = rtol * resnorm + atol
    ItData = ItStatsPTC(resnorm)
    #
    # Preallocate a vector for the step
    #
    step = copy(x)
    #
    # If the initial iterate satisfies the termination criteria, tell me.
    #
    toosoon = (resnorm <= tol)
    #
    # The main loop stops on convergence or too many iterations.
    #
    delta = delta0
    while resnorm > tol && itc < maxit
        #
        # Evaluate and factor the Jacobian; update x, F(x), and delta.
        #
        (x, delta, FS, resnorm) = PTCUpdate(FPS, FS, x, ItRules, step, resnorm, delta)
        #
        # Keep the books
        #
        updateStats!(ItData, resnorm)
        itc += 1
        ~keepsolhist || (@views solhist[:, itc+1] .= x)
    end
    (idid, errcode) = PTCOK(resnorm, tol, toosoon, ItRules, printerr)
    itout = CloseIteration(x, FS, ItData, idid, errcode, keepsolhist, solhist)
    return (itout)
end


text:
# helper function for 'RT-locking' and 'stretching' of GC-program
function stretched_program(n::Float64, par::GasChromatographySimulator.Parameters)
	# stretch the temperature program in 'par' by a factor 'n'
	if isa(par, Array)==true
		error("Select an element of the array of GC-system parameters.")
	else
		new_tsteps = n.*par.prog.time_steps
		new_T_itp = GasChromatographySimulator.temperature_interpolation(new_tsteps, par.prog.temp_steps, par.prog.gf, par.col.L)
		new_pin_itp = GasChromatographySimulator.pressure_interpolation(new_tsteps, par.prog.pin_steps)
		new_pout_itp = GasChromatographySimulator.pressure_interpolation(new_tsteps, par.prog.pout_steps)
		new_prog = GasChromatographySimulator.Program(new_tsteps, par.prog.temp_steps, par.prog.pin_steps, par.prog.pout_steps, par.prog.gf, par.prog.a_gf, new_T_itp, new_pin_itp, new_pout_itp)
		new_par = GasChromatographySimulator.Parameters(par.col, new_prog, par.sub, par.opt)
		return new_par
	end
end

function initial_n(n::Float64, tR_lock::Float64, ii::Int, par::GasChromatographySimulator.Parameters)
    par_n = stretched_program(n, par)
    sol_n = GasChromatographySimulator.solve_system_multithreads(par_n)
    tR_n = sol_n[ii].u[end][1]
    if n>1
        while tR_n-tR_lock<0 && n<130.0
            n = n*2.0
            par_n = stretched_program(n, par)
            sol_n = GasChromatographySimulator.solve_system_multithreads(par_n)
            tR_n = sol_n[ii].u[end][1]
        end
    elseif n<1
        while tR_n-tR_lock>0 && n>0.01
            n = n*0.5
            par_n = stretched_program(n, par)
            sol_n = GasChromatographySimulator.solve_system_multithreads(par_n)
            tR_n = sol_n[ii].u[end][1]
        end
    end
    if n>130.0
        error("The choosen retention time for locking is to big.")
    elseif n<0.01
        error("The choosen retention time for locking is to small.")
    else
        return n
    end
end

function RT_locking(par::GasChromatographySimulator.Parameters, tR_lock::Float64, tR_reltol::Float64, solute_RT::String; opt_itp="linear")
	# estimate the factor 'n' for the temperature program to achieve the retention time 'tR_lock' for 'solute_RT' with the GC-system defined by 'par'
	if isa(par, Array)==true
		error("Select an element of the array of GC-systems.")
	elseif tR_reltol<par.opt.reltol
		error("The tolerance for retention time locking tR_reltol is to small. Use tR_reltol >= par.opt.reltol ($(par.opt.reltol)).")
	else
		# find 'solute_RT' in the substances of 'par'
		name = Array{String}(undef, length(par.sub))
		for i=1:length(par.sub)
			name[i] = par.sub[i].name
		end
		ii = findfirst(name.==solute_RT)
		# calculate the retention time for the original (un-stretched) program
		sol₀ = GasChromatographySimulator.solving_odesystem_r(par.col, par.prog, par.sub[ii], par.opt)
		tR₀ = sol₀.u[end][1]
		# start value for the factor 'n'
		if tR₀-tR_lock<0
			n₁ = initial_n(2.0, tR_lock, ii, par)
		else
			n₁ = initial_n(0.5, tR_lock, ii, par)
		end
		# using a recursive function to estimate 'n'
		n = recur_RT_locking(n₁, [1.0], [tR₀], par, tR_lock, tR_reltol, ii; opt_itp=opt_itp)
	end
	return n
end

function recur_RT_locking(n::Float64, n_vec::Array{Float64,1}, tR_vec::Array{Float64,1}, par::GasChromatographySimulator.Parameters, tR_lock::Float64, tR_reltol::Float64, ii::Int64; opt_itp="linear")
	# recursive function to find the factor 'n' for the temperature program to achieve the retention time 'tR_lock' for solute index 'ii' with the GC-system defined by 'par'
	# calculate the retention time with the input guess 'n'
	par₁ = stretched_program(n, par)
	if par.opt.odesys==true
		sol₁ = GasChromatographySimulator.solving_odesystem_r(par₁.col, par₁.prog, par₁.sub[ii], par₁.opt)
		tR₁ = sol₁.u[end][1]
	else
		sol₁ = GasChromatographySimulator.solving_migration(par₁.col, par₁.prog, par₁.sub[ii], par₁.opt)
		tR₁ = sol₁.u[end]
	end
	if abs(tR₁-tR_lock)/tR_lock<tR_reltol
		# if retention time is less than 'tR_tol' from 'tR_lock' we found the factor 'n'
		return n
	else
		# estimate a new factor 'new_n' by linear interpolation of the factors 'n_vec' + 'n' over the corresponding retention times
		new_n_vec = sort([n_vec; n])
		new_tR_vec = sort([tR_vec; tR₁])
		if opt_itp=="spline"
			# Dierckx.jl
			if length(new_tR_vec)<4
				k = length(new_tR_vec)-1
			else
				k = 3
			end
			itp = Spline1D(new_tR_vec, new_n_vec, k=k)
		else # opt_itp=="linear"
			# Interpolations.jl
			itp = LinearInterpolation(sort([tR_vec; tR₁]), sort([n_vec; n]))
		end
		new_n = itp(tR_lock)
		#println("new_n=$(new_n), tR₁=$(tR₁)")
		# use the new factor 'new_n' and call the recursive function again
		return recur_RT_locking(new_n, new_n_vec, new_tR_vec, par, tR_lock, tR_reltol, ii; opt_itp=opt_itp)
	end
end


text:
function _exportdensematrix!(
    A::GBVecOrMat{T};
    desc = nothing
) where {T}
desc = _handledescriptor(desc)
    nrows = Ref{libgb.GrB_Index}(size(A,1))
    ncols = Ref{libgb.GrB_Index}(size(A,2))
    Csize = Ref{libgb.GrB_Index}(length(A) * sizeof(T))
    values = Ref{Ptr{Cvoid}}(Ptr{T}())
    isuniform = Ref{Bool}(false)
    libgb.GxB_Matrix_export_FullC(
        Ref(A.p),
        Ref(toGBType(T).p),
        nrows,
        ncols,
        values,
        Csize,
        isuniform,
        desc
    )
    A.p = C_NULL
    finalize(A)
    return nrows[], ncols[], values[]

end
function _exportdensematrix(A::GBMatrix; desc = nothing)
    return _exportdensematrix!(copy(A); desc)
end
function Base.Matrix(A::GBMatrix{T}) where {T}
    if gbget(A, SPARSITY_STATUS) != GBDENSE
        X = similar(A)
        X[:, :] = zero(T)
        A = eadd(X, A)
    end
    nrows, ncols, values = _exportdensematrix(A)
    C = Matrix{T}(undef, nrows, ncols)
    unsafe_copyto!(pointer(C), Ptr{T}(values), length(C))
    ccall(:jl_free, Cvoid, (Ptr{T},), values)
    return C
end

function Vector(v::GBVector{T}) where {T}
    if gbget(v, SPARSITY_STATUS) != GBDENSE
        X = similar(v)
        X[:] = zero(T)
        v = eadd(X, v)
    end
    nrows, _, vals = _exportdensematrix!(copy(v))
    v = Vector{T}(undef, nrows)
    unsafe_copyto!(pointer(v), Ptr{T}(vals), length(v))
    ccall(:jl_free, Cvoid, (Ptr{T},), vals)
    return v
end

function _exportcscmatrix!(
    A::GBMatrix{T};
    desc = nothing
    ) where {T}
    desc = _handledescriptor(desc)
    nrows = Ref{libgb.GrB_Index}(size(A, 1))
    ncols = Ref{libgb.GrB_Index}(size(A, 2))
    t = Ref{libgb.GrB_Type}(toGBType(T).p)
    colptr = Ref{Ptr{libgb.GrB_Index}}()
    rowidx = Ref{Ptr{libgb.GrB_Index}}()
    values = Ref{Ptr{Cvoid}}(Ptr{T}())
    colptrsize = Ref{libgb.GrB_Index}()
    rowidxsize = Ref{libgb.GrB_Index}()
    Axsize = Ref{libgb.GrB_Index}()
    isuniform = Ref{Bool}(false)
    isjumbled = C_NULL

    libgb.GxB_Matrix_export_CSC(
        Ref(A.p),
        t,
        nrows,
        ncols,
        colptr,
        rowidx,
        values,
        colptrsize,
        rowidxsize,
        Axsize,
        isuniform,
        isjumbled,
        desc
    )
    A.p = C_NULL
    finalize(A)
    nrows = nrows[]
    ncols = ncols[]
    colptr = colptr[]
    rowidx = rowidx[]
    colptrsize = colptrsize[]
    rowidxsize = rowidxsize[]
    Axsize = Axsize[]
    val = values[]
    return nrows, ncols, colptr, rowidx, colptrsize, rowidxsize, val, Axsize

end

function _exportcscmatrix(A::GBMatrix; desc = nothing)
    return _exportcscmatrix!(copy(A); desc)
end

function SparseArrays.SparseMatrixCSC(A::GBMatrix{T}; desc = nothing) where {T}
    nrows, ncols, colptr, rowidx, colptrsize, rowidxsize, val, valsize = _exportcscmatrix(A; desc)
    outvalues = Vector{T}(undef, valsize ÷ sizeof(T))
    col = Vector{libgb.GrB_Index}(undef, Int(colptrsize ÷ sizeof(libgb.GrB_Index)))
    row = Vector{libgb.GrB_Index}(undef, Int(rowidxsize ÷ sizeof(libgb.GrB_Index)))
    unsafe_copyto!(pointer(outvalues), Ptr{T}(val), length(outvalues))
    unsafe_copyto!(pointer(col), Ptr{libgb.GrB_Index}(colptr), length(col))
    unsafe_copyto!(pointer(row), Ptr{libgb.GrB_Index}(rowidx), length(row))
    ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), colptr)
    ccall(:jl_free, Cvoid, (Ptr{libgb.GrB_Index},), rowidx)
    ccall(:jl_free, Cvoid, (Ptr{T},), val)
    return SparseArrays.SparseMatrixCSC(nrows, ncols, col .+= 1, row .+= 1, outvalues)
end

function _exportdensevec!(
    v::GBVector{T};
    desc = nothing
    ) where {T}
    desc = _handledescriptor(desc)
    n = Ref{libgb.GrB_Index}(size(v,1))
    vsize = Ref{libgb.GrB_Index}(length(v) * sizeof(T))
    values = Ref{Ptr{Cvoid}}(Ptr{T}())
    isuniform = Ref{Bool}(false)
    libgb.GxB_Vector_export_Full(
        Ref(v.p),
        Ref(toGBType(T).p),
        n,
        values,
        vsize,
        isuniform,
        desc
    )
    return n[], values[]
end

function _exportdensevec(v::GBVector; desc = nothing)
    return _exportdensevec!(copy(v); desc)
end