#include <math.h>

#include "ElementCheckFunction.h"
//#include "CheckElement_Unit.h"
#include "ProcExt.h"
#include "ListConst.h"
//#include "CheckParam.h"
//#include "bkerip_ext.h"
#include "sigma_beam.h"

//------------------------------------------------------------------------------
// Вычисление коэффициента fi_b для центрально-сжатых элементов
double Get_fi_b_centr(int type_sect, COMMON_SECT *com_sect, double M_, double N_, double len_design_YZ,
                double length_ext, double length_int,
                GEOM_PARAM_SECT *param_sect, double Ry, double Ey) {

       double tw, hw, bf, tf, hb;
       double alfa_176, psi_, fi_1;
       double Force[7] = {0, 0, 0, 0, 0, 0, 0};
       double sigma_comm[4];
       double len_ef;

       double fi_b;

       if (type_sect != DVUTAVR)  {
         fi_b = 1;
         return fi_b;
       }

       tw = com_sect->dvutavr.b;  // толщина стенки
       hw = com_sect->dvutavr.h;  // высота стенки
       bf = com_sect->dvutavr.b1;  // ширина полки
       tf = com_sect->dvutavr.h1;  // толщина полки  (внешняя и внутренняя полка имеют одинаковую толщину)
       hb = hw + tf; // расстояние между осями поясов

       //--------------------------------------
       // относительный эксцентриситет
       Force[0] = N_;
       Force[2] = M_;

       sigma_beam_comm(param_sect, Force, sigma_comm);
       if (sigma_comm[0] < sigma_comm[2]) {
         if (length_ext == 0)
           len_ef = len_design_YZ;
         else
           len_ef = length_ext;
       }
       else {
         if (length_int == 0)
           len_ef = len_design_YZ;
         else
           len_ef = length_int;
       }
        // 2. Вычисление fi_b     // стр. 85. Приложение 7*
       if (com_sect->dvutavr.n_group < 0) { // сварной двутавр
          alfa_176 = Round((8.0*pow(len_ef*tf/(hb*bf),2.0))*(1 + 0.5*hb*tw*tw*tw/(bf*tf*tf*tf)), 2);   // ф.176
       }
       else {   // Прокатной двутавр
          alfa_176 = Round(1.54*param_sect->Itor*pow(len_ef/param_sect->Hsect, 2.0)/param_sect->Iyy, 2);   // ф.176
       }

       psi_ = Round(Fun_psi (alfa_176), 2);             // табл. 77
       fi_1 = Round((psi_*param_sect->Iyy/param_sect->Izz)*pow(param_sect->Hsect/len_ef, 2.0)*Ey/Ry, 2); // ф. 174
       if (fi_1 <=0.85)
          fi_b = 0.85;
       else
          fi_b = Round((MIN(1.0, 0.68 + 0.21*fi_1)), 2);
       return fi_b;
}
//------------------------------------------------------------------------------
// Вычисление коэффициента fi_b для изгибаемых элементов
// по приложению Ж
double Get_fi_b_bend_comp(int type_sect, COMMON_SECT *com_sect, double width, double thick, double M_, double N_, double len_design_YZ,
				double length_ext, double length_int,
				GEOM_PARAM_SECT *param_sect, double Ry, double Ey,
				char type_constr, char num_constr, char type_load, char index_poyas0, FI_B_PARAM *Fi_b_param, int index_comp) {

	 COMMON_SECT com_sect_comp;
	 double Af2;
	 double fi_b;

	 memcpy((char*)&com_sect_comp, (char*)com_sect, sizeof(COMMON_SECT));

	 switch (type_sect) {
	  case DVUTAVR:
	  Af2 = com_sect->dvutavr.h2*com_sect->dvutavr.b2 + width*thick;
	  com_sect_comp.dvutavr.b2 = width;
	  com_sect_comp.dvutavr.h2 = Af2/width;
	 break;
	 case TAVR:
	  if (index_comp==0) {
		com_sect_comp.dvutavr.h2 = thick;
		com_sect_comp.dvutavr.b2 = width;
	  }
	  else {
		com_sect_comp.dvutavr.h2 = com_sect_comp.dvutavr.h1;
		com_sect_comp.dvutavr.b2 = com_sect_comp.dvutavr.b1;
		com_sect_comp.dvutavr.h1 = thick;
		com_sect_comp.dvutavr.b1 = width;
	  }
	 break;
	 }
	 fi_b = Get_fi_b_bend(DVUTAVR, &(com_sect_comp), M_, N_,
						  len_design_YZ, length_ext, length_int,
						  param_sect, Ry, Ey,
						  type_constr, num_constr, type_load, index_poyas0, Fi_b_param);

	 return fi_b;
}
//------------------------------------------------------------------------------
// Вычисление коэффициента fi_b для изгибаемых элементов
// по приложению Ж
double Get_fi_b_bend(int type_sect, COMMON_SECT *com_sect, double M_, double N_, double len_design_YZ,
				double length_ext, double length_int,
				GEOM_PARAM_SECT *param_sect, double Ry, double Ey,
				char type_constr, char num_constr, char type_load, char index_poyas0, FI_B_PARAM *Fi_b_param) {

   // type_constr;  // тип закрепления 0 - пролет, 1 - консоль
   // num_constr;   // количество закреплений в пролете 0 - без закреплений, 1 - два в середине 2 - одно в середине
   // type_load;    // тип нагрузки 0 - сосредоточенная в середине, 1 - равномерно распределенная; 2 - сосредоточенная на четверти 2 - любая;
   // index_poyas;  // 0 - нагрузка приложена к сжатому поясу, 1 - к растянутому

     double tw, hw, bf, tf, hb;
     double alfa_, psi, psi_1, fi_1;
     double Force[7] = {0, 0, 0, 0, 0, 0, 0};
     double sigma_comm[4];
     double len_ef;
     int index_comp;  // 0 - сжат наружный пояс, 1 - внутренний
     char index_poyas;

     double fi_b;

     index_poyas = MIN(1, index_poyas0);

     if (type_sect != DVUTAVR)  {
       fi_b = 1;
       return fi_b;
     }
     Force[0] = N_;
     Force[2] = M_;

	 //if (Round(M_, 0) == 268374667) {
     //  fi_b = 1;
     //}

     sigma_beam_comm(param_sect, Force, sigma_comm);
     if (sigma_comm[0] < sigma_comm[2]) {
	   index_comp = 0;
       if (length_ext == 0)
         len_ef = len_design_YZ;
       else
         len_ef = length_ext;
     }
     else {
       index_comp = 1;
       if (length_int == 0)
         len_ef = len_design_YZ;
       else
         len_ef = length_int;
     }

     tw = com_sect->dvutavr.b;  // толщина стенки
     hw = com_sect->dvutavr.h;  // высота стенки
	 bf = com_sect->dvutavr.b1;  // ширина полки
	 if (com_sect->dvutavr.h1 == com_sect->dvutavr.h2 && com_sect->dvutavr.b1 == com_sect->dvutavr.b2 ) { // двутавр с двумя осями симметрии
	   tf = com_sect->dvutavr.h1;  // толщина полки

	   hb = hw + tf; // расстояние между осями поясов

		// 2. Вычисление fi_b     // СП стр. 141. Приложение Ж
	   if (com_sect->dvutavr.n_group < 0) { // сварной двутавр
		  alfa_ = Round((8.0*pow(len_ef*tf/(hb*bf),2.0))*(1 + 0.5*hb*tw*tw*tw/(bf*tf*tf*tf)), 2);   // Ж.5
	   }
	   else {   // Прокатной двутавр
		  alfa_ = Round(1.54*param_sect->Itor*pow(len_ef/param_sect->Hsect, 2.0)/param_sect->Iyy, 2);   // Ж.4
	   }

	   // Назначение psi
	   if (type_constr==0) {  // пролет
		 switch (num_constr) {
		   case 0:
			 if (type_load==0) { // Сосредоточенная
			   if (index_poyas==0) { // К сжатому
				 if (alfa_<=40)
				   psi = 1.75 + 0.09*alfa_;
				 else
				   psi = 3.3 +0.053*alfa_ - 4.5*1e-5*alfa_*alfa_;
			   }
			   else {  // К растянутому
				 if (alfa_<=40)
				   psi = 5.05 + 0.09*alfa_;
				 else
				   psi = 6.6 +0.053*alfa_ - 4.5*1e-5*alfa_*alfa_;
			   }
			 }
			 else {  // Равномерно распределенная и другие
			   if (index_poyas==0) { // К сжатому
				 if (alfa_<=40)
				   psi = 1.6 + 0.08*alfa_;
				 else
				   psi = 3.15 + 0.04*alfa_ - 2.7*1e-5*alfa_*alfa_;
			   }
			   else {  // К растянутому
				 if (alfa_<=40)
				   psi = 3.8 + 0.08*alfa_;
				 else
				   psi = 5.35 + 0.04*alfa_ - 2.7*1e-5*alfa_*alfa_;
			   }
			 }
		   break;
		   case 1:   // количество закреплений сжатого пояса два и более
			  if (alfa_<=40)
				 psi = 2.25 + 0.07*alfa_;
			  else
				 psi = 3.6 + 0.04*alfa_ - 3.5*1e-5*alfa_*alfa_;
		   break;
		   case 2:   // количество закреплений сжатого пояса одно в середине
			  if (alfa_<=40)
				 psi_1 = 2.25 + 0.07*alfa_;
			  else
				 psi_1 = 3.6 + 0.04*alfa_ - 3.5*1e-5*alfa_*alfa_;
			 if (type_load==0) { // Сосредоточенная в середине
			   psi = 1.75*psi_1;
			 }
			 if (type_load==1) { // равномерно распределенная
			   if (index_poyas==0)  // К сжатому
				  psi = 1.14*psi_1;
			   else
				  psi = 1.30*psi_1;
			 }
			 if (type_load==2) { // сосредоточенная в четверти
			   if (index_poyas==0)  // К сжатому
				  psi = 1.14*psi_1;
			   else
				  psi = 1.60*psi_1;
			 }
		   break;
		   default:
			  fi_b = 1;
			  return fi_b;
		 }
	   }
	   else {  // Консоль
		  if (type_load==0) { // Сосредоточенная
			   if (index_poyas==0) { // К сжатому
				 if (alfa_<=28)
				   psi = 6.2 + 0.08*alfa_;
				 else
				   psi = 7.0 + 0.05*alfa_;
			   }
			   else {  // К растянутому
				 if (alfa_<=28)
				   psi = 1.00 + 0.16*alfa_;
				 else
				   psi = 4.0 + 0.05*alfa_;
			   }
		  }
		  else {   // Равсномерно распределенная
			 psi = 1.42*sqrt(alfa_);
		  }

	   }
	   double h_fi_b;
	   if (com_sect->dvutavr.n_group < 0) // Сварной двутавр
		 h_fi_b = param_sect->Hsect - com_sect->dvutavr.h1;
	   else
		 h_fi_b = param_sect->Hsect;

	   fi_1 = Round((psi*param_sect->Iyy/param_sect->Izz)*pow(h_fi_b/len_ef, 2.0)*Ey/Ry, 2); // ф. Ж.3

	   if (fi_1 <=0.85)
		  fi_b = fi_1;
	   else
		  fi_b = Round((MIN(1.0, 0.68 + 0.21*fi_1)), 2);
	   Fi_b_param->psi = psi;
	   Fi_b_param->fi_1 = fi_1;
	   Fi_b_param->h_fi_b = h_fi_b;
	   Fi_b_param->l_ef = len_ef;
	 }

	 else {  // двутавр с одной осью симметрии
	   double beta, eta, delta, mu;
	   double B_, C_, D_, C09, C10, k_, psi_a;
	   double n0_, fi_1, fi_2;
	   double tf1, tf2, b1, b2, h1, h2;
	   double I1, I2;

	   tf1 = com_sect->dvutavr.h1;  // толщина полки внешнего пояса
	   tf2 = com_sect->dvutavr.h2;  // толщина полки внутреннего пояса
	   b1 = com_sect->dvutavr.b1;
	   b2 = com_sect->dvutavr.b2;
	   tf = (tf1 + tf2)*0.5;
	   hb = hw + tf; // расстояние между осями поясов

	   alfa_ = Round((8.0*pow(len_ef*tf/(hb*bf),2.0))*(1 + 0.5*hb*tw*tw*tw/(bf*tf*tf*tf)), 2);   // Ж.5

	   // Вычисление B, C и D
	   if (tf1*b1 > tf2*b2) {
		 h1 = param_sect->Zct - tf1*0.5;
		 h2 = hb - h1;
		 I1 = tf1*b1*b1*b1/12 + tf1*b1*(hw + tf1)*(hw + tf1)/4;
		 I2 = tf2*b2*b2*b2/12 + tf2*b2*(hw + tf2)*(hw + tf2)/4;

		 n0_ = I1/(I1 + I2);
		 beta = (2*n0_ - 1)*(0.47 - 0.035*(b1/hb)*(1 + b1/hb - 0.072*(pow(b1/hb, 2.0))));
		 eta = (1 - n0_)*(9.87*n0_ + 0.385*(param_sect->Itor/I2)*pow(len_ef/hb, 2.0));
		 delta = n0_ + 0.734*beta;
		 mu = n0_ + 1.145*beta;

		 if (index_comp==0) { // сжат внешний и он более развит
			if (index_poyas == 0) { // нагрузка приложена к сжатому
			  if (type_load==0)  B_= delta - 1;
			  else B_ = mu - 1;
			}
			else {
			  if (type_load==0)  B_= delta;
			  else B_ = mu;
			}
		 }
		 else {  //сжат внутренний и он менее развит
			if (index_poyas == 0) { // нагрузка приложена к сжатому
			  if (type_load==0)  B_= -delta;
			  else B_ = -mu;
			}
			else {
			  if (type_load==0)  B_= 1 - delta;
			  else B_ = 1 - mu;
			}
		 }
	   }
	   else {
		 h2 = param_sect->Zct - tf1*0.5;
		 h1 = hb - h2;
		 I2 = tf1*b1*b1*b1/12 + tf1*b1*(hw + tf1)*(hw + tf1)/4;
		 I1 = tf2*b2*b2*b2/12 + tf2*b2*(hw + tf2)*(hw + tf2)/4;
		 n0_ = I1/(I1 + I2);

		 beta = (2*n0_ - 1)*(0.47 - 0.035*(b1/hb)*(1 + b1/hb - 0.072*(pow(b1/hb, 2.0))));
		 eta = (1 - n0_)*(9.87*n0_ + 0.385*(param_sect->Itor/I2)*pow(len_ef/hb, 2.0));
		 delta = n0_ + 0.734*beta;
		 mu = n0_ + 1.145*beta;

		 if (index_comp==0) { // сжат внешний и он менее развит
			if (index_poyas == 0) { // нагрузка приложена к сжатому
			  if (type_load==0)  B_= -delta;
			  else B_ = -mu;
			}
			else {
			  if (type_load==0)  B_= 1 - delta;
			  else B_ = 1 - mu;
			}
		 }
		 else {  //сжат внутренний и он более развит
			if (index_poyas == 0) { // нагрузка приложена к сжатому
			  if (type_load==0)  B_= delta - 1;
			  else B_ = mu - 1;
			}
			else {
			  if (type_load==0)  B_= delta;
			  else B_ = mu;
			}
		 }
	   }
	   if (type_load==0) {
		 C09 = 0.33*eta;
		 C10 = 0.0826*alfa_;
		 D_ = 3.265;
	   }
	   else {
		 C09 = 0.481*eta;
		 C10 = 0.1202*alfa_;
		 D_ = 2.247;
	   }
	   if (n0_ <= 0.9)
		 C_ = C09;
	   else {
		 if (n0_ < 1.0) {
		   C_ = C09 + (C10 - C09)*(n0_ - 0.9)/0.1;
		 }
		 else
		   C_ = C10;
	   }

	   psi_a = (B_ + sqrt(B_*B_ + C_))*D_;

	   fi_1 = Round((psi_a*param_sect->Iyy/param_sect->Izz)*(2*hb*h1/(len_ef*len_ef))*Ey/Ry, 2); // ф. Ж.6
	   fi_2 = Round((psi_a*param_sect->Iyy/param_sect->Izz)*(2*hb*h2/(len_ef*len_ef))*Ey/Ry, 2); // ф. Ж.7
	   if ((tf1 > tf2 && index_comp==1) || (tf1 < tf2 && index_comp==0)) {  // Сжатый пояс менее развит
		  if (n0_ >0.7 && len_ef/bf >= 5 && len_ef/bf <= 25) {
			k_ = 1.025 - 0.015*len_ef/bf;
			fi_2 *= k_;
			fi_2 = MIN(0.95, fi_2);
		  }
		  if (fi_2 <= 0.85)
			fi_b = fi_2;
		  else
			fi_b = MIN(1.0, 0.68 + 0.21*fi_2);
	   }
	   else {  // сжатый пояс более развит
		  if (fi_2 <= 0.85)
			fi_b = MIN(1.0, fi_1);
		  else
			fi_b = MIN(1.0, fi_1*(0.21+0.68*(n0_/fi_1 + (1 - n0_)/fi_2)));
	   }
	   Fi_b_param->psi_a = psi_a;
	   Fi_b_param->fi_1 = fi_1;
	   Fi_b_param->fi_2 = fi_2;
	   Fi_b_param->l_ef = len_ef;
	   Fi_b_param->h_fi_b = hb;

	 }
	 if (fi_b==0) {
	   fi_b = 1;
	 }
	 Fi_b_param->fi_b = fi_b;
	 return fi_b;
}
