//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop

#include "ChoiceKolonna_Alg.h"
#include "ListConst.h"
#include "message.h"
#include "ProcExt.h"
#include "Data_Gen.h"
#include "StandartProfil.h"
#include "Dvutavr_Choice_Alg.h"
#include "Get_param_sect_func.h"



//---------------------------------------------------------------------------
#pragma package(smart_init)
extern MATER MaterProp[MAX_MATERIAL];
extern int   Material_Number;
//---------------------------------------------------------------------------
// Выполнить подбор сечения для колонн постоянного сечения
int ChoiсeKolonna_Alg(PARAM_CHOSEN_SECT *Param_chosen_sect, SECT_DVUTAVR *sect,
					  double lx_ef, double ly_ef, double *mat_stress, bool flag_simm) {
	 GEOM_PARAM_SECT param_sect;
	 double kR;
	 double sig_ext, sig_int;
     double lambda_w, lambda_x, lambda_;
	 double fi_e, h_opt, hw, tw, tf, bf, mult_ix, ix;
	 int i;
	 double gamma_c;
	 double e_, Bx_, By_;
     double dzeta, omega, Const0, Const1, lambda_f;
     double lambda_y, Rel_u, ksi_, dzeta_sk, u_, uef, Cx_, Delta, Delta_e;
     double Ar_, Rel_A, eta, omega_ef, k_;
	 double m_exc, lam, Add_koef, Up_lambda_w, ss;
	 int n_loop;
	 double eta_work[5];
	 TStandartProfil StandartProfil;
	 MATER_PARAM mater_param;
	 double Ry, Rs, Rsw, Ey, Ryw;
	 double thick_range[2];
	 int rc;
	 SECT_DVUTAVR sect_dubl;
	 double Mz, Nx;


	 gamma_c = Param_chosen_sect->gamma_c;
	 double min_hsten = Param_chosen_sect->min_hsten;
	 double max_bpol = Param_chosen_sect->max_bpol;
	 double min_bpol = Param_chosen_sect->min_bpol;
	 double max_tpol = Param_chosen_sect->max_tpol;
	 double min_tpol = Param_chosen_sect->min_tpol;
	 double max_tsten = Param_chosen_sect->max_tsten;
	 double min_tsten = Param_chosen_sect->min_tsten;

	 double round_h = Param_chosen_sect->round_h;
	 double round_width = Param_chosen_sect->round_width;
	 double ind_round_thick = Param_chosen_sect->round_thick;

	 /*   тест для гибкой стенки
	 mu_x = 2.0;
	 mu_y = 0.8;
	 l_kol =8600;
	 Nx = 259294;
	 Mz = 71761000;
	 gamma_c =0.95;
	 */
	 /* тест для устойчивой стенки
	 mu_x = 2.0;
	 mu_y = 0.4615;
	 l_kol =7800;
	 Nx = 930000;
	 Mz = 319000000;
	 gamma_c =1.0;
	 Ry = 210;
	 */

	 Mz = ABS(mat_stress[1]);
	 Nx = ABS(mat_stress[0]);
	 Nx = MAX(1, Nx);

	 double Mx=2*Mz/3;
	 // 1. Предварительный расчет
	 // Эксцентриситет
	 e_ = Mz/Nx;
	 if (e_> 1e5) {
		return ERR_ALG2_BEND;
	 }

	 char stroka[240];
	 #ifdef PRINT_DEBUG
	 InfMessage("Alg_Gepard-----------------");
	 sprintf(stroka, "lx_ef=%10.2f, ly_ef=%10.2f, e_=%10.2f", lx_ef, ly_ef, e_);
	 InfMessage((AnsiString)stroka);
	 #endif
	 n_loop =0;
	 Ey = MaterProp[Material_Number].E;

loop0:
	// Расчетное сопротивление растяжению, сжатию и изгибу по пределу текучески
	 double t_max = MAX(sect->h1, sect->h2);
	//  Получить таблицу свойств стали по марке стали
	 rc = Get_Mater_param(MaterProp[Material_Number].title, &mater_param);
	 if (rc>0) return rc;
	// Проверить границы толщин и получить диапазон толщин
	 rc = Get_thick_range(MaterProp[Material_Number].title, t_max, thick_range);
	 if (rc==1 || rc==2 || rc==11) {  // толщины не попадают в диапазон толщин, но ничего не предпринимать
		 ; //return 100;
	 }
	 Ry = Steel_param_Ry(MaterProp[Material_Number].title, thick_range[0], &rc, true);

	 // Тест !!!
	 //Ry=210;
	 //Ey=206000;

	 // Вычисление характерных параметров - перенесено за метку, т. Ry может уточняться (пока не сделано)
	 kR = sqrt(Ey/Ry);

	 Bx_ = Nx*Ey/(Ry*gamma_c*Ry*lx_ef*lx_ef);

	 dzeta = lx_ef/(2*ly_ef);

	 omega = 1.3*e_* kR/lx_ef;

	 By_ = Nx*Ey/(Ry*gamma_c*Ry*ly_ef*ly_ef);

	 Const0 = 5.0*(1+0.26*dzeta*dzeta)*(1+0.26*dzeta*dzeta)*sqrt(Ry/Ey)/(1+1.8*omega);
	 Const1 = (1+0.26*dzeta*dzeta)*(1+0.26*dzeta*dzeta)*sqrt(Ry/Ey)/(1+2.9*omega);

	 if (By_ > Const0) {
	   lambda_ = 1.0;
	   By_ = Const0;
	 }
	 else {
	   if (By_ <= Const0 && By_ >= Const1) lambda_ = 2.0;
	   if (By_ < Const1) lambda_ = 3.0;
	 }
	 if (Param_chosen_sect->flag_gib)  // гибкая стенка
	   lambda_y = 1.5*lambda_;
	 else
	   lambda_y = lambda_;

	 // Конец вычисления характерных параметров
	 // Начало итераций
	 t_max = thick_range[0];
	 bool flag_loop = true;
loop:
	// Расчетное сопротивление растяжению, сжатию и изгибу по пределу текучески
	//  Получить таблицу свойств стали по марке стали
	 rc = Get_Mater_param(MaterProp[Material_Number].title, &mater_param);
	 if (rc>0) return rc;
	// Проверить границы толщин и получить диапазон толщин
	 rc = Get_thick_range(MaterProp[Material_Number].title, t_max, thick_range);
	 if (rc==1 || rc==2 || rc==11) {  // толщины не попадают в диапазон толщин
     	 Ry = Steel_param_Ry(MaterProp[Material_Number].title, t_max, &rc, true);
		 return 100;
	 }
	 Ry = Steel_param_Ry(MaterProp[Material_Number].title, t_max, &rc, true);

	 // Тест !!!
	 // Ry=210;
	 // Ey=206000;

	 // Гибкость поясов
	 lam = MAX(0.8, lambda_);
	 lam = MIN(4.0, lambda_);
	 lambda_f = 2.0*(0.36 + 0.1*lam)*kR;

	 if (Param_chosen_sect->flag_gib) { // гибкая стенка
		Add_koef = 1.7;
		Up_lambda_w = 4.6*kR;
        if (lambda_<2.0)
          lambda_w = Add_koef*(1.3 +0.15*lambda_*lambda_)*kR;
		else
		  lambda_w = Add_koef*(1.2 +0.35*lambda_)*kR;
        lambda_w = MIN(Up_lambda_w, lambda_w);
		//lambda_y = 1.5*lambda_;
		Rel_u = 1.7*kR/lambda_w;
		ksi_ = 3.0*Rel_u*(1 - Rel_u + Rel_u*Rel_u/3);
        ss = sqrt(2.0*ksi_*lambda_f*pow(lx_ef*lambda_y/(ly_ef*lambda_),2)/(9.0*lambda_w) + 1);
        dzeta_sk = sqrt((3.0 * lambda_w/(ksi_*lambda_f))*
                   (ss - 1));
        u_ = lambda_f*dzeta_sk*dzeta_sk/(2*lambda_w);
		uef = u_ * Rel_u;
	 }
	 else {
		if (lambda_<2.0)
          lambda_w = (1.3 +0.15*lambda_*lambda_)*kR;
        else
          //lambda_w = (1.2 +0.35*lambda_)*kR;
          lambda_w = (MIN(2.3, (1.2 +0.35*lambda_)))*kR;
		//lambda_y = lambda_;
        ksi_ = 1;
        ss = sqrt(2.0*ksi_*lambda_f*pow(lx_ef*lambda_y/(ly_ef*lambda_),2)/(9.0*lambda_w) + 1);
        dzeta_sk = sqrt((3.0 * lambda_w/(ksi_*lambda_f))*
				   (ss - 1));
		u_ = lambda_f*dzeta_sk*dzeta_sk/(2*lambda_w);
		uef = u_;
	 }

	 Cx_ = u_*(3.0+ksi_*u_)*lambda_w/(12.0*(1+uef)*(1+uef));
	 Delta = Bx_*Cx_;
	 omega  = sqrt(3.0*(1+uef)/(3 + ksi_*u_))*e_*kR/lx_ef;
	 Delta_e = (1+3*pow(omega,0.75)/(4*pow(Delta,0.3333)))*Delta;
		  // Уточнение Delta_e
	 Delta_e = (1+omega*Delta_e/(pow(Delta_e+0.04, 1.5)))*Delta;
	 Delta_e = MAX(0.03, Delta_e);

	 lambda_x = 1/(sqrt(Delta_e+0.04));
	 Ar_ = Nx *(1+omega*Delta_e/(pow(Delta_e+0.04, 1.5)))*
		   (1 + 0.04/Delta_e)/(Ry*gamma_c);

	 #ifdef PRINT_DEBUG
	 sprintf(stroka, "lambda_f=%10.2f, lambda_y=%10.2f", lambda_f, lambda_y);
	 InfMessage((AnsiString)stroka);
	 sprintf(stroka, "lambda_w=%10.2f, u_=%10.2f, uef=%10.2f", lambda_w, u_, uef);
	 InfMessage((AnsiString)stroka);
	 sprintf(stroka, "Delta_e=%10.2f, lambda_x=%10.2f, Ar_=%10.2f", Delta_e, lambda_x, Ar_);
	 InfMessage((AnsiString)stroka);
	 #endif
	 if ((ABS(lambda_x - lambda_))/lambda_ >0.05 && n_loop<5) {  // Повторение расчета
		lambda_y = lambda_y*lambda_x/lambda_;
		lambda_ = (lambda_ + 2*lambda_x)/3;
		//if (n_loop==0) {      // !!!!!!!!!!!!!!!!!! - отладка
		//  lambda_y = 1.957;
		//  lambda_ = 1.638;
		//}
		n_loop++;
		goto loop;
	 }
	 else {
		lambda_ = lambda_x;
	 }
 //---------------------------------------------
 // Вычисление требуемой площади сечения
 // вычисление относительного эксцентриситета
	 m_exc = omega*lambda_;
	 m_exc = Round(m_exc, 2);

	 Rel_A = 1/(2*u_);
	 eta = Fun_eta(Rel_A, lambda_, m_exc, eta_work);

	 omega_ef = eta * omega;

	 if (Delta>0.1) k_=0.16*(2+Delta)/(0.3+Delta);
	 if (Delta<=0.1 && Delta>=0.03) k_=0.84;
	 if (Delta<0.03) k_=0.51 + 11*Delta;
	 if (Delta>=0.03) {
		 lambda_x = 1/sqrt((1+k_*pow(omega_ef,0.8))*(0.035+Delta));
	 }
	 else {
		 lambda_x = 1.67/(sqrt((1+k_*pow(omega_ef,0.8))*sqrt(Delta)));
	 }
	 fi_e = Delta*lambda_x*lambda_x;

	 Ar_ = Nx/(fi_e*Ry*gamma_c);        // требуемая площадь сечения

	 #ifdef PRINT_DEBUG
	 sprintf(stroka, "2: lambda_x=%10.2f, Ar_=%10.2f", lambda_x, Ar_);
	 InfMessage((AnsiString)stroka);
	 #endif

	 // Компоновка сечения
	 mult_ix = sqrt((3+uef)/(12*(1+uef)));

	 ix = lx_ef*sqrt(Ry/Ey)/lambda_x;   // радиус инерции

	 h_opt = ix/mult_ix;                // высота сечения

	 bf = h_opt/dzeta_sk;               // ширина полки

	 tf = bf/lambda_f;                  // толщина полки
	 tw = (h_opt - 2*tf)/lambda_w;      // толщина стенки

	 if (Param_chosen_sect->flag_minimax) {   // Учитывать ограничения на размеры
	   h_opt = MAX(h_opt - 2*tf, min_hsten) + 2*tf;
	 // Повторное назначение размеров сечения
	   bf = h_opt/dzeta_sk;               // ширина полки
	   bf = MAX(bf, min_bpol);               // ширина полки
	   bf = MIN(bf, max_bpol);               // ширина полки

	   tf = bf/lambda_f;
	   tf = MAX(tf, min_tpol);                  // толщина полки
	   tf = MIN(tf, max_tpol);                  // толщина полки

	   tw = (h_opt - 2*tf)/lambda_w;      // толщина стенки
	   tw = MAX(tw, min_tsten);
	   tw = MIN(tw, max_tsten);
	 }

	 memcpy(&(sect_dubl), sect, sizeof(SECT_DVUTAVR));

	 //--------------------------------------------------
	 // Чтобы вычислить геометрические параметры сечения
	 hw = Round_step(h_opt - 2*tf, Param_chosen_sect->round_h);
	 sect->h = hw;
	 sect->b = tw;
	 sect->b1 = bf;
	 sect->h1 = tf;
	 sect->b2 = bf;
	 sect->h2 = tf;
	 // вычисление геометрических характеристик сечения
	 Get_param_sect(sect, &(param_sect));
	 double Izz0 = param_sect.Izz;

	 if (sect->n_group < 0 ) {   // сварное сечение
	   tf = Round_row(tf, ind_round_thick);
	   if (ind_round_thick == 2) {
		 for (i=0; i<4; i++) {
		   double bef0 = Get_max_bpol(tf, ind_round_thick);
		   if (bf > bef0) {
			 sect->b1 = bef0;
			 sect->h1 = tf;
			 sect->b2 = bef0;
			 sect->h2 = tf;
			 Get_param_sect(sect, &(param_sect));
			 if (param_sect.Izz < Izz0*0.99)
			   tf = Round_row(tf+1, ind_round_thick);
			 else
			   break;
		   }
		 }
	   }
	   bf = Round_row_bpol(tf, bf, Param_chosen_sect->round_width, ind_round_thick);
	   tw = Round_row(tw, ind_round_thick);
	   hw = Round_step(h_opt - 2*tf, Param_chosen_sect->round_h);

	   sect->h = hw;
	   sect->b = tw;
	   sect->b1 = bf;
	   sect->h1 = tf;
	   sect->b2 = bf;
	   sect->h2 = tf;
	 }
	 if (sect->n_group >= 0 ) {   // стандартное сечение
		//int rc = StandartProfil.Choice_Profil_param(sect, Ar_, ix, 0.0, h_opt, tw, tf, bf, true);
		int rc = StandartProfil.Choice_Profil_param(sect, Ar_, 0.0, 0.0, h_opt, tw, tf, bf, true);
		if (rc>0)
		  return rc;
	 }
	 t_max = MAX(sect->h1, sect->h2);
	 if (flag_loop) {
		 flag_loop = false;
		 goto loop;
	 }
	 // вычисление геометрических характеристик сечения
	 Get_param_sect(sect, &(param_sect));
	 // напряжения во внешнем поясе
	 sig_ext = -param_sect.Zct*(Mz/param_sect.Izz) + Nx/param_sect.Ar;
	 // напряжения во внутреннем поясе
	 sig_int = (param_sect.Hsect - param_sect.Zct)*(Mz/param_sect.Izz) + Nx/param_sect.Ar;

	 return 0;

}
//----------------------------------------------------------------------------
// Функция вычисления коэффициента eta для двутавра  (табл. 73-5)
double Fun_eta(double A_, double lambda_, double m_, double *eta_work) {
	   double eta_025, eta_050, eta_100;
	   if (lambda_<=5.0) {
          if (m_<=5.0) {
             eta_025 = (1.45 - 0.05*m_) - 0.01*(5-m_)*lambda_;
             eta_050 = (1.75 - 0.10*m_) - 0.02*(5-m_)*lambda_;
             eta_100 = (1.90 - 0.10*m_) - 0.02*(6-m_)*lambda_;
          }
		  else {
             if (m_<=20) {
                eta_025 = 1.2;
                eta_050 = 1.25;
                eta_100 = 1.4 - 0.02*lambda_;
             }
			 else {
                eta_025 = 1.2;
                eta_050 = 1.25;
                eta_100 = 1.4 - 0.02*lambda_;
             }
          }
       }
       else {
          eta_025 = 1.2;
          eta_050 = 1.25;
          eta_100 = 1.3;
       }
       if (A_<=0.25) return eta_025;
       if (A_>=1.0)  return eta_100;
       eta_work[0] = 0;
       eta_work[1] = 0;
       eta_work[2] = 0;
       eta_work[3] = 0;
       if (A_<=0.5) {
         eta_work[0] = eta_025;
         eta_work[1] = eta_050;
         eta_work[2] = 0.25;
         eta_work[3] = 0.5;
		 return eta_025 + (A_ - 0.25)*(eta_050 - eta_025)/0.25;
	   }
       else {
         eta_work[0] = eta_050;
         eta_work[1] = eta_100;
         eta_work[2] = 0.5;
         eta_work[3] = 1.0;
         return eta_050 + (A_ - 0.50)*(eta_100 - eta_050)/0.50;  // if (A_<1.0)
       }
}


