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

#include <math.h>

#pragma hdrstop

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

static double Tabl_74[21][25]= {SNIP_Tabl_74};
static double Line_74[25]= {SNIP_Line_74};
static double Column_74[21]= {SNIP_Column_74};
static double Tabl_18[6][8]= {SP_Tabl_18};
static double Line_18[8]= {SP_Line_18};
static double Column_18[6]= {SP_Column_18};
static double delta_cr[7]={0.8, 1., 2.0, 4.0, 6.0, 10.0, 30.0};
static double c_cr[7]={30.0, 31.5, 33.3, 34.6, 34.8, 35.1, 35.5};
static double Tabl_14[7][10]={SP_Tabl_14};
static double Line_14[10] ={SP_Line_14};
static double Column_14[7] ={SP_Column_14};
static double Tabl_15[6][8]={SP_Tabl_15};
static double Line_15[8] ={SP_Line_15};
static double Column_15[6] ={SP_Column_15};
//------------------------------------------------------------
// Вычисление силы Qfic
double GetQfic(double Ey, double Ry, double N_, double fi) {
   return (7.15*1e-6*(2330 - Ey/Ry)*((ABS(N_))/fi));
}
//---------------------------------------------------------------------------
//  Функция определения fi_e аппроксимацией по матрице  (табл. 74 СНиП)
double Fun_fi_e(double lambda_r0, double m_ef0, double fi_) {
   double fi_e, di, dj, g1, g2;
   double m_ef, lambda_r;
   int i1, j1;
       lambda_r = MIN(lambda_r0, 14.0);
       m_ef = MIN(m_ef0, 20.0);

	   fi_e = App_matrix_fun(Tabl_74[0], 21, 25, Column_74, Line_74, m_ef, lambda_r);

	   fi_e = MIN(fi_e, fi_);

       return fi_e;
}
//---------------------------------------------------------------------
//  Функция определения alfa аппроксимацией по матрице  (табл. 18 СП)
double Fun_alfa_18(double lambda_w0, double tau_R0) {
   double alfa, di, dj, g1, g2;
   double tau_R, lambda_w;
   int i1, j1;
       lambda_w = MIN(lambda_w0, 5.5);
	   tau_R = MIN(tau_R0, 0.9);

       alfa = App_matrix_fun(Tabl_18[0], 6, 8, Column_18, Line_18, lambda_w, tau_R);
       /*
			for (i1=0; i1<6; i1++)
			  if (Column_18[i1]>tau_R) break;
            i1--;
			for (j1=0; j1<8; j1++)
              if (Line_18[j1]>lambda_w) break;
			j1--;
            i1=MIN(6,i1);
            j1=MIN(8,j1);
			if (i1==-1) {
              di = 0.0;
              i1=0;
            }
			else
              di = (tau_R - Column_18[i1])/(Column_18[i1+1] - Column_18[i1]);
            if (j1==-1) {
              j1 = 0;
			  dj = 0;
            }
            else
              dj = (lambda_w - Line_18[j1])/(Line_18[j1+1] - Line_18[j1]);
            g1 = Tabl_18[i1][j1]*(1 - di) + Tabl_18[i1+1][j1]*di;
            g2 = Tabl_18[i1][j1+1]*(1 - di) + Tabl_18[i1+1][j1+1]*di;
			alfa = g1*(1 - dj) + g2*dj;
			*/
            return alfa;
}
//----------------------------------------------------------------------------
// Функция вычисления коэффициента alfa
double Fun_alfa(double mx) {
    if (mx <= 1.0)
       return 0.7;
    else
       return(0.65 + 0.05*mx);   //     if (mx <= 5.0)
}
//----------------------------------------------------------------------------
// Функция вычисления коэффициента beta
double Fun_beta(double lambda_y, double lambda_c, double fi_y, double fi_c) {
    if (lambda_y <= lambda_c)
	   return 1.0;
    else
       return (sqrt(fi_c/fi_y));
}
//----------------------------------------------------------------------------
// Функция вычисления коэффициента psi, табл. 77 СНиП
double Fun_psi(double alfa_) {
    if (alfa_ <= 40.0)
       return 2.25 + 0.07*alfa_;
	else
       return 3.6 + 0.04*alfa_ -3.5e-5*alfa_*alfa_;
}
//----------------------------------------------------------------------------
// Функция вычисления коэффициента eta для двутавра  по оси стенки (табл. 73-п.5 или Д.2 - п.5)
double Fun_eta_dvutavr(double A_, double lambda_, double m_, double *eta_work) {
       // A_ - отношение площади поясов к площади стенки
       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)
	   }
}
//----------------------------------------------------------------------------
// Функция вычисления коэффициента eta для двутавра  по оси стенки (табл. 73-п.5, п.10 или Д.2 - п.5, п.10)
double Fun_eta_dvutavr5_10(double Af1, double Af2, double Aw, double lambda_, double m_, double *eta_work) {
	   // A_ - отношение площади поясов к площади стенки
	   double eta, eta5, eta10;
	   double A_, af;
	   double Af_max, Af_min;
	   double eta_025, eta_050, eta_100, eta_200;
	   Af_max = MAX(Af1, Af2);
	   Af_min = MIN(Af1, Af2);
	   A_ = Af_max/Aw;
	   // Вычисление для двутавра с одинаковыми поясами, рывными максимальному
	   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) eta5 = eta_025;
	   if (A_>=1.0)  eta5 = 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;
		 eta5 = 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;
		 eta5 = eta_050 + (A_ - 0.50)*(eta_100 - eta_050)/0.50;  // if (A_<1.0)
	   }
	   if (Af1 != Af2) {  // п.10
		 if (lambda_<=5.0) {
		   if (m_<=5.0) {
			 eta_050 = 1.4;
			 eta_100 = 1.6 - 0.01*(5 - m_)*lambda_;
			 eta_200 = 1.8 - 0.02*(5 - m_)*lambda_;
		   }
		   else {
			 eta_050 = 1.4;
			 eta_100 = 1.6;
			 eta_200 = 1.8;
		  }
		 }
		 else {
		   if (m_<=5.0) {
			 eta_050 = 1.4;
			 eta_100 = 1.35 + 0.05*m_;
			 eta_200 = 1.3 + 0.1*m_;
		   }
		   else {
			 eta_050 = 1.4;
			 eta_100 = 1.6;
			 eta_200 = 1.8;
		   }
		 }
		 if (A_<=0.5) eta10 = eta_050;
		 if (A_>=2.0) eta10 = eta_200;

		 if (A_<=1.0) {
		   eta10 = eta_050 + (A_ - 0.5)*(eta_100 - eta_050)/0.5;
		 }
		 else {
		   eta10 = eta_100 + (A_ - 1.00)*(eta_200 - eta_100)/1.00;  // if (A_<2.0)
		 }
		 af = Af_min/Af_max;
		 eta = eta5 + (1.0 - af)*(eta10 - eta5)/0.5;
	   }
	   else {
		 eta = eta5;
	   }
	   return eta;
}
//----------------------------------------------------------------------------
// Функция вычисления коэффициента eta для двутавра по оси перпендикулярной стенки (табл. 73-п.8 или Д.2 - п.8)
double Fun_eta_dvutavr8(double A_, double lambda_, double m_, double *eta_work) {
	   // A_ - отношение площади стенки к площади поясов
	   double eta_025, eta_050, eta_100;
	   if (lambda_<=5.0) {
		  if (m_<=5.0) {
			 eta_025 = (0.75 + 0.05*m_) + 0.01*(5-m_)*lambda_;
			 eta_050 = (0.50 + 0.10*m_) + 0.02*(5-m_)*lambda_;
			 eta_100 = (0.25 + 0.15*m_) + 0.03*(5-m_)*lambda_;
		  }
		  else {
			 if (m_<=20) {
				eta_025 = 1.0;
				eta_050 = 1.0;
				eta_100 = 1.0;
             }
             else {
                eta_025 = 1.0;
                eta_050 = 1.0;
				eta_100 = 1.0;
             }
          }
       }
	   else {
          eta_025 = 1.0;
          eta_050 = 1.0;
          eta_100 = 1.0;
       }
       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)
       }
}
//----------------------------------------------------------------------------
// Функция вычисления коэффициента eta для тавра  (табл. 73-11)
double Fun_eta_tavr_w(double A_, double lambda_, double m_, double *eta_work) {
       // A_ - отношение площади поясов к площади стенки
	   double eta_050, eta_100, eta_150, eta_200;
	   if (lambda_<=5.0) {
          if (m_<=5.0) {
             eta_050 = (1.45 + 0.04*m_);
             eta_100 = (1.80 + 0.12*m_);
             eta_150 = (2.00 + 0.25*m_) + 0.1*lambda_;
             eta_200 = (3.00 + 0.25*m_) + 0.1*lambda_;
          }
          else {
             eta_050 = 1.65;
             eta_100 = 2.4;
             eta_150 = 100;
             eta_200 = 100;
          }
       }
       else {
         if (m_<=5.0) {
			 eta_050 = (1.45 + 0.04*m_);
             eta_100 = (1.80 + 0.12*m_);
             eta_150 = 100;
             eta_200 = 100;
         }
         else {
             eta_050 = 1.65;
             eta_100 = 2.4;
             eta_150 = 100;
			 eta_200 = 100;
         }
       }
       if (A_<=0.5) return eta_050;
	   if (A_>=2.0) return eta_200;
       eta_work[0] = 0;
	   eta_work[1] = 0;
       eta_work[2] = 0;
       eta_work[3] = 0;
       if (A_<=1.0) {
         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.5)*(eta_100 - eta_050)/0.5;
       }
       else {
         if (A_<1.5) {
		   eta_work[0] = eta_100;
		   eta_work[1] = eta_150;
		   eta_work[2] = 1.0;
		   eta_work[3] = 1.5;
		   return eta_100 + (A_ - 1.0)*(eta_150 - eta_100)/0.50;  // if (A_<1.5)
		 }
		 else {
		   eta_work[0] = eta_150;
		   eta_work[1] = eta_200;
		   eta_work[2] = 1.5;
		   eta_work[3] = 2.0;
		   return eta_150 + (A_ - 1.5)*(eta_200 - eta_150)/0.50;  // if (A_<2.0)
		 }
	   }
}
//----------------------------------------------------------------------------
// Функция вычисления коэффициента eta (табл. Д.2 (73-5 СНиП))
double Fun_eta_sect(int type_sect, COMMON_SECT *com_sect, double lambda_, double m_, double *eta_work, int index_Y) {
   double eta_, Af, Af1, Af2, Aw;
   switch (type_sect) {
	 case DVUTAVR:
	  Af1 = com_sect->dvutavr.h1*com_sect->dvutavr.b1;
	  Af2 = com_sect->dvutavr.h2*com_sect->dvutavr.b2;
	  Af = (Af1 + Af2)/2;

	  Aw = com_sect->dvutavr.h*com_sect->dvutavr.b;
	  if (index_Y==0) {
		//eta_ = Fun_eta_dvutavr(Af/Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5
		eta_ = Fun_eta_dvutavr5_10(Af1, Af2, Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5-10
	  }
	  else
		eta_ = Fun_eta_dvutavr8(Aw/(2*Af), lambda_, m_, eta_work);   // табл. Д2, тип 8
	 break;
	 case RECT_TUBE:
	  Af = com_sect->rect.b*com_sect->rect.b1;
	  Aw = com_sect->rect.h*com_sect->rect.h1;
	  if (index_Y==0) {
		eta_ = Fun_eta_rect_tube(Af, 2*Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5-10
	  }
	  else
		eta_ = Fun_eta_rect_tube(2*Af, Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5-10
	 break;
	 case SQUARE_TUBE:
	  Af = com_sect->square.side*com_sect->square.thickness;
	  Aw = Af;
	  if (index_Y==0) {
		eta_ = Fun_eta_rect_tube(Af, 2*Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5-10
	  }
	  else
		eta_ = Fun_eta_rect_tube(2*Af, Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5-10
	 break;
	 case CIRCULAR_TUBE:
	  eta_ = Fun_eta_circ_tube(lambda_, m_);   // табл. 73-4
	  eta_work[0] = 0;
	  eta_work[1] = 0;
	  eta_work[2] = 0;
	  eta_work[3] = 0;
	 break;
	 //case TAVR:
	 case TAVR_UP:
	 case TAVR_DOWN:
	  Af = (com_sect->dvutavr.h1*com_sect->dvutavr.b1 + com_sect->dvutavr.h2*com_sect->dvutavr.b2)/2;
	  Aw = com_sect->dvutavr.h*com_sect->dvutavr.b;
	  eta_ = Fun_eta_tavr_w(Af/Aw, lambda_, m_, eta_work);   // табл. 73-11
	 break;
	 case RECTANCULAR:
	   eta_ = 1.0;
	 break;
	 default:
	  //int i = Application->MessageBox(("Fun_eta_sect:Сечение типа "+IntToStr(type_sect)+
	  //                                 " не обрабатывается, принимается eta=1").c_str(),
	  //            " ", MB_OK | MB_ICONWARNING);

	  eta_ = 1.0;
   }
   return eta_;
}
//----------------------------------------------------------------------------
// Функция вычисления коэффициента eta (табл. Д.2 (73-5 СНиП))  с учетом пластины
double Fun_eta_sect_plate(int type_sect, COMMON_SECT *com_sect, double lambda_, double m_, double *eta_work, int index_Y,
						 double plate_w, double plate_h, int index_comp) {
   double eta_, Af, Af1, Af2, Aw;
   switch (type_sect) {
	 case DVUTAVR:
	  Af1 = com_sect->dvutavr.h1*com_sect->dvutavr.b1;
	  Af2 = com_sect->dvutavr.h2*com_sect->dvutavr.b2;
	  if (index_comp>=0) {
		  Af2 += plate_w*plate_h;
	  }
	  Af = (Af1 + Af2)/2;

	  Aw = com_sect->dvutavr.h*com_sect->dvutavr.b;
	  if (index_Y==0) {
		//eta_ = Fun_eta_dvutavr(Af/Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5
		eta_ = Fun_eta_dvutavr5_10(Af1, Af2, Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5-10
	  }
	  else
		eta_ = Fun_eta_dvutavr8(Aw/(2*Af), lambda_, m_, eta_work);   // табл. Д2, тип 8
	 break;
	 case RECT_TUBE:
	  Af = com_sect->rect.b*com_sect->rect.b1;
	  Aw = com_sect->rect.h*com_sect->rect.h1;
	  if (index_Y==0) {
		eta_ = Fun_eta_rect_tube(Af, 2*Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5-10
	  }
	  else
		eta_ = Fun_eta_rect_tube(2*Af, Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5-10
	 break;
	 case SQUARE_TUBE:
	  Af = com_sect->square.side*com_sect->square.thickness;
	  Aw = Af;
	  if (index_Y==0) {
		eta_ = Fun_eta_rect_tube(Af, 2*Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5-10
	  }
	  else
		eta_ = Fun_eta_rect_tube(2*Af, Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5-10
	 break;
	 case CIRCULAR_TUBE:
	  eta_ = Fun_eta_circ_tube(lambda_, m_);   // табл. 73-4
	  eta_work[0] = 0;
	  eta_work[1] = 0;
	  eta_work[2] = 0;
	  eta_work[3] = 0;
	 break;
	 case TAVR:
	 case TAVR_UP:
	 case TAVR_DOWN:
	  Aw = com_sect->tavr.h*com_sect->tavr.b;
	  if (index_comp>=0) {
		 Af1 = com_sect->tavr.h1*com_sect->tavr.b1;
		 Af2 = plate_w*plate_h;
		 eta_ = Fun_eta_dvutavr5_10(Af1, Af2, Aw, lambda_, m_, eta_work);   // табл. Д2, тип 5-10
	  }
	  else {
		Af = com_sect->tavr.h1*com_sect->tavr.b1;
		eta_ = Fun_eta_tavr_w(Af/Aw, lambda_, m_, eta_work);   // табл. 73-11
	  }
	 break;
	 case RECTANCULAR:
	   eta_ = 1.0;
	 break;
	 default:
	  //int i = Application->MessageBox(("Fun_eta_sect:Сечение типа "+IntToStr(type_sect)+
	  //                                 " не обрабатывается, принимается eta=1").c_str(),
	  //            " ", MB_OK | MB_ICONWARNING);

	  eta_ = 1.0;
   }
   return eta_;
}
//----------------------------------------------------------------------------
// Функция вычисления коэффициента eta для трубы  (табл. Д2)
double Fun_eta_circ_tube(double lambda_, double m_) {
	double eta_;
	if (lambda_<=5.0) {
	   if (m_<=5.0) {
		  eta_ = (1.35 - 0.05*m_) - 0.01*(5-m_)*lambda_;
	   }
	   else
		  eta_ = 1.1;
	}
	else {
	   eta_ = 1.1;
	}
	return eta_;
}
//----------------------------------------------------------------------------
// Функция вычисления коэффициента eta для прямоугольной трубы  (Д2)
double Fun_eta_rect_tube(double Af, double Aw, double lambda_, double m_, double *eta_work) {
	   // A_ - отношение площади поясов к площади стенки
	   double eta, eta5, eta10;
	   double A_, af;
	   double eta_025, eta_050, eta_100, eta_200;
	   A_ = Af/Aw;
	   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) eta5 = eta_025;
	   if (A_>=1.0)  eta5 = 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;
		 eta5 = 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;
		 eta5 = eta_050 + (A_ - 0.50)*(eta_100 - eta_050)/0.50;  // if (A_<1.0)
	   }
	   return eta5;
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------
// Функция вычисления коэффициента fi по формулам 8, 9, п.7.1.3 СП
double Fun_fi_713(int type_sect, double lambda_r, int index_Y) {
	double fi_, delta;
	double alfa, beta;
	double max_lambda_r;

	if (lambda_r >= 0.4) {
	  switch (type_sect) {
	  case CIRCULAR_TUBE:
	  case SQUARE_TUBE:
	  case RECT_TUBE:
       alfa = 0.03;
       beta = 0.06;
       max_lambda_r = 3.8;
      break;
	  case RECTANCULAR:
	   alfa = 0.04;
	   beta = 0.09;
	   max_lambda_r = 4.4;
	  break;
	  case DVUTAVR:
	   alfa = 0.04;
	   if (index_Y==0)
		  beta = 0.09;  // В плоскости стенки
	   else
		  beta = 0.14;  // В плоскости поясов
	   max_lambda_r = 4.4;
	  break;
      default:
       alfa = 0.04;
       beta = 0.14;
       max_lambda_r = 5.8;
	  }
	  delta = 9.87*(1 - alfa + beta*lambda_r) + lambda_r*lambda_r;
      fi_ = 0.5*(delta - sqrt(delta*delta - 39.48*lambda_r*lambda_r))/(lambda_r*lambda_r);
	  if (lambda_r > max_lambda_r)
		 fi_ = MIN(7.6/(lambda_r*lambda_r), fi_);
    }
	else
      fi_ = 1;
    return fi_;
}
//---------------------------------------------------------------------
// Функция вычисления условной предельной гибкости стенки при центральном сжатии по табл. 9 СП
double Fun_lambda_uw_cent(int type_sect, double lambda_r) {
    double lambda_uw;

      switch (type_sect) {
      case CIRCULAR_TUBE:
      case SQUARE_TUBE:
      case RECT_TUBE:
       if (lambda_r <= 1.0) {
         lambda_uw = 1.2;
       }
       else {
         lambda_uw = MIN(1.6, 1.0 + 0.2*lambda_r);
       }
      break;
      case DVUTAVR:
       if (lambda_r <= 2.0) {
         lambda_uw = 1.3 + 0.15*lambda_r*lambda_r;
       }
	   else {
         lambda_uw = MIN(2.3, 1.2 + 0.35*lambda_r);
       }
      break;
	  case SHVELLER:
       if (lambda_r <= 0.8) {
         lambda_uw = 1.0;
	   }
       else {
         lambda_uw = MIN(1.6, 0.85 + 0.19*lambda_r);
       }
      break;
      default:
		lambda_uw = (0.4 + 0.07*lambda_r)*(1+0.25);
		#ifdef EMBRAC
		int i = Application->MessageBox((L"Fun_lambda_uw_cent:Сечение типа "+IntToStr(type_sect)+
										 " не обрабатывается, принимается lambda_uw="+FloatToStr(lambda_uw)).c_str(),
				  L" ", MB_OK | MB_ICONWARNING);
		#else
		int i = Application->MessageBox(("Fun_lambda_uw_cent:Сечение типа "+IntToStr(type_sect)+
										 " не обрабатывается, принимается lambda_uw="+FloatToStr(lambda_uw)).c_str(),
				  " ", MB_OK | MB_ICONWARNING);
		#endif
	  }
      return lambda_uw;
}
//--------------------------------------------------------------------------------------
// Функция вычисления предельной гибкости поясов при центральном сжатии по табл. 10 СП
double Fun_lambda_uf_cent(int type_sect, double lambda_r) {
    double lambda_uf, lambda_r0;

      lambda_r0 = MAX(0.8, lambda_r);
      lambda_r0 = MIN(4.0, lambda_r);

	  switch (type_sect) {
	  case CIRCULAR_TUBE:
      case SQUARE_TUBE:
	  case RECT_TUBE:   // по табл. 9
	   if (lambda_r <= 1.0) {
		 lambda_uf = 1.2;
       }
       else {
         lambda_uf = MIN(1.6, 1.0 + 0.2*lambda_r);
       }
	  break;
      case DVUTAVR:
        lambda_uf = 0.36 + 0.1*lambda_r0;
      break;
      default:
        lambda_uf = (0.43 + 0.08*lambda_r0);
      }
      return lambda_uf;
}
//--------------------------------------------------------------------------------------
// Функция вычисления предельной гибкости стенки при mx = 20 по п. 8.5.8 СП
double Fun_lambda_uw20(double M_, double Q_, double Ryf,  double Ryw, double Ey, double Rsw, double hef,
					   double t_, double Af, double gamma_c) {
   double r_, alfa_f, alfa, tau;
   double ti, krit, lambda_w;
   int i;

   ti = t_;
   r_ = MAX(1.0, Ryf/Ryw);
   for (i=0; i<10; i++) {
	 lambda_w = MIN(5.5, (hef/ti)*sqrt(Ryw/Ey));        // Табл. 18
     if ((ABS(M_)) <1e-1)
	   break;
     alfa_f = Af/(ti*hef);
     tau = (ABS(Q_))/(ti*hef);
	 alfa = Fun_alfa_18(lambda_w, tau/Rsw);
     krit = (ABS(M_))/(Ryf*gamma_c*hef*hef*ti*(r_*alfa_f + alfa));   // из ф. 86
     if (krit < 1.0 && lambda_w == 5.5)
       break;
     ti = ti*krit;
   }
   return lambda_w;
}
//--------------------------------------------------------------------------------------
// Функция вычисления предельной гибкости поясов при mx = 20 по п. 8.5.18 или 8.5.19 СП
double Fun_lambda_uf20(int index_105, double Ryf,  double sigma_c, double lambda_uw) {
   if (index_105==1) { // п. 8.5.19 - пластика допускается для балок 2 и 3 классов
      return (0.17 + 0.06*lambda_uw);
   }
   else {  // п. 8.5.19 - пластика НЕ допускается для балок 1 класса
	  return (0.5*sqrt(Ryf/(ABS(sigma_c))));
   }
}
//--------------------------------------------------------------------------------------
// Функция вычисления предельной гибкости поясов при mx = 20 по п. 8.5.18 или 8.5.19 СП
double Fun_rect_lambda_uf20(int index_105, double Ryf,  double sigma_c, double lambda_uw) {
   if (index_105==1) { // п. 8.5.19 - пластика допускается для балок 2 и 3 классов
	  return (0.675 + 0.15*lambda_uw);
   }
   else {  // п. 8.5.19 - пластика НЕ допускается для балок 1 класса
	  return (1.5*sqrt(Ryf/(ABS(sigma_c))));
   }
}
//--------------------------------------------------------------------------------------
// Функция вычисления предельной гибкости стенки по табл. 22 СП - при внецентренном сжатии
double Fun_lambda_uw_cent_bend(int type_sect, double lambda_x, double m_exc, double sig1, double sig2, double tau,
							   double delta_fi, double Ry, double gamma_c, double lambda_uw20)
 {
	  double lambda_uw, lambda_uw0, m_exc5, lambda_uw1;
	  double c_cr_17[6] = {C_CR_17};
	  double alfa_cr_17[6] = {ALFA_CR_17};
	  double alfa, c_cr, beta;
	  double lambda_uw_05, lambda_uw_1;
	  double alfa_1;
	  double sig1_plus;

	  lambda_uw0 = Fun_lambda_uw_cent(type_sect, lambda_x);

      if (sig1!=0.0)
        alfa = (sig1 - sig2)/sig1;
      else
        alfa = 1;  
      switch (type_sect) {
      case CIRCULAR_TUBE:
      case SQUARE_TUBE:
      case RECT_TUBE:
        if (lambda_x < 2.0)
          lambda_uw1 = lambda_uw = 1.3 + 0.15*lambda_x*lambda_x;
        else
		  lambda_uw1 = lambda_uw = MIN (1.2 + 0.35*lambda_x, 3.1);
        if (m_exc < 1) {
            lambda_uw0 = Fun_lambda_uw_cent(type_sect, lambda_x);
            lambda_uw = lin_interpol(0.0, 1.0, lambda_uw0, lambda_uw1, m_exc);
        }
        if (m_exc > 10 && m_exc < 20) {
            lambda_uw = lin_interpol(10.0, 20.0, lambda_uw1, lambda_uw20, m_exc);
		}
        if (m_exc >= 20) {
            lambda_uw = lambda_uw20;
        } 
	  break;
      case DVUTAVR:
        if (delta_fi > 0) {   // delta_fi = c_*fi_y - fi_e, т.е. c_*fi_y >= fi_e
          if (lambda_x < 2.0)
			lambda_uw1 = 1.3 + 0.15*lambda_x*lambda_x;
          else
            lambda_uw1 = MIN (1.2 + 0.35*lambda_x, 3.1);
		  lambda_uw = lambda_uw1;     // m_xc >= 1.0 && m_exc <= 10
          if (m_exc < 1) {
            lambda_uw0 = Fun_lambda_uw_cent(type_sect, lambda_x);
            lambda_uw = lin_interpol(0.0, 1.0, lambda_uw0, lambda_uw1, m_exc);
          }
          if (m_exc > 10 && m_exc < 20) {
            lambda_uw = lin_interpol(10.0, 20.0, lambda_uw1, lambda_uw20, m_exc);
          }
          if (m_exc >= 20) {
            lambda_uw = lambda_uw20;
          }
		  //InfMessage("c*fi_y > fi_e: delta_fi=" + FloatToStrF(delta_fi, ffFixed, 15, 4)+"; lambda_uw="+ FloatToStrF(lambda_uw, ffFixed, 15, 2));
        }
        else {  // c_*fi_y <= fi_e
          sig1_plus = ABS(sig1);
          if (alfa >= 1) {        // 1<=alfa<=2 - alfa не может быть больше двух, т.к. при sig1 = -sig2 alfa = 2
            c_cr = l_app(alfa_cr_17, c_cr_17, alfa, 6);
            beta = 0.15*c_cr*tau/sig1_plus;
            lambda_uw = MIN(1.42*sqrt(c_cr*Ry*gamma_c/(sig1_plus*(2 - alfa + sqrt(alfa*alfa + 4*beta*beta)))),
							0.7 + 2.4*alfa);
          }
          else {
             lambda_uw_05 = MIN(Fun_lambda_uw_cent(type_sect, lambda_x),
							  Fun_lambda_uw_cent_bend(type_sect, lambda_x, m_exc, sig1, sig2, tau,
                              1.0, Ry, gamma_c, lambda_uw20));
             if (alfa >=0.5) {
               alfa_1 = 1.0;
               c_cr = l_app(alfa_cr_17, c_cr_17, alfa_1 , 6);
			   beta = 0.15*c_cr*tau/sig1_plus;
               lambda_uw_1 = MIN(1.42*sqrt(c_cr*Ry*gamma_c/(sig1_plus*(2 - alfa_1 + sqrt(alfa_1*alfa_1 + 4*beta*beta)))),
                            0.7 + 2.4*alfa);
			   lambda_uw = lambda_uw_05 + (lambda_uw_1 - lambda_uw_05)*(1 - alfa)/0.5;
            }
            else {
              lambda_uw = lambda_uw_05;
            }
          }
          //InfMessage("c*fi_y <= fi_e: delta_fi=" + FloatToStrF(delta_fi, ffFixed, 15, 4)+"; lambda_uw="+ FloatToStrF(lambda_uw, ffFixed, 15, 2));
          //InfMessage("sig1=" + FloatToStrF(sig1, ffFixed, 15, 1)+"; sig2="+ FloatToStrF(sig2, ffFixed, 15, 1)+"; alfa="+FloatToStrF(alfa, ffFixed, 15, 2));
          //InfMessage("c_cr=" + FloatToStrF(c_cr, ffFixed, 15, 2)+"; beta="+ FloatToStrF(beta, ffFixed, 15, 2));
        }
      break;
      default:
         lambda_uw = 1;
      }
      return lambda_uw;
}
//--------------------------------------------------------------------------------------
// Функция вычисления предельной гибкости поясов по табл. 23 СП - при внецентренном сжатии
double Fun_lambda_uf_cent_bend(int type_sect, double lambda_cent, double lambda_r, double lambda_uf20, double m_exc, double Ry) {
	double lambda_uf, lambda_c0, lambda_r0, lambda_ufc, m_exc5, lambda_uf5;

	  lambda_c0 = MAX(0.8, lambda_cent);
	  lambda_c0 = MIN(4.0, lambda_cent);

	  lambda_ufc = Fun_lambda_uf_cent(type_sect, lambda_c0);

	  lambda_r0 = MAX(0.8, lambda_r);
	  lambda_r0 = MIN(4.0, lambda_r);

      m_exc5 = MIN(m_exc, 5.0);
      switch (type_sect) {
	  case CIRCULAR_TUBE:
      case SQUARE_TUBE:
	  case RECT_TUBE:
		lambda_uf5 = lambda_ufc - 0.01*(5.3 + 1.3*lambda_r0)*m_exc5;
      break;
      case DVUTAVR:
		lambda_uf5 = lambda_ufc - 0.01*(1.5 + 0.7*lambda_r0)*m_exc5;
      break;
      default:
		lambda_uf5 = lambda_ufc - 0.01*(1.5 + 0.7*lambda_r0)*m_exc5;
      }
      if (m_exc >= 20)
        return lambda_uf20;
      if (m_exc > 5) {
        lambda_uf = lambda_uf5 + (m_exc - 5)*(lambda_uf20 - lambda_uf5)/15;
        return lambda_uf;
      }
      else
        return lambda_uf5;
}
//---------------------------------------------------------------------
// Получить значение коэффициентов n, cx, cy по табл Е.1  СП 16.13330.2011
void Get_n_cx_cy(int index_axis, int type_sect, COMMON_SECT *com_sect, double *n_, double *cx_, double *cy_) {

	  double af_aw_1[4] = {AF_AW_1};    // Для  симметричного сечения
	  double af_aw_2[3] = {AF_AW_2};    // Для  неимметричного двутавра
	  double af, aw, cx_1, cx_2, cy_1, cy_2, n_1, n_2;
	  double Af1, Af2, Af_min, Af_max;
	  double Cx_1[4] = {CX_dvutavr_1};  // Для  симметричного двутавра
	  double Cx_2[3] = {CX_dvutavr_2};  // Для  несимметричного двутавра

	  double Cx_rect[4] = {CX_rect};  // Для  прямоугольного сечения
	  double Cy_rect[4] = {CY_rect};  // Для  прямоугольного сечения

	  switch (type_sect) {
	  case SQUARE_TUBE:
		*cx_ = 1.12;
		*cy_ = 1.12;
		*n_ = 1.5;

	  break;
	  case RECT_TUBE:
		if (index_axis == 0) { // в плоскости Y-Y
		  Af1 = com_sect->rect.b*com_sect->rect.b1;
		  aw = 2*com_sect->rect.h*com_sect->rect.h1;

		  cx_1 = l_app(af_aw_1, Cx_rect, Af1/aw, 4);
		  cy_1 = l_app(af_aw_1, Cy_rect, Af1/aw, 4);
		  n_1 = 1.5;
		  *cx_ = Round(cx_1, 3);
		  *cy_ = Round(cy_1, 3);
		  *n_ = Round(n_1, 3);
		}
		else { // в плоскости X-X
		  Af1 = com_sect->rect.h*com_sect->rect.h1;
		  aw = 2*com_sect->rect.b*com_sect->rect.b1;

		  cx_1 = l_app(af_aw_1, Cx_rect, Af1/aw, 4);
		  cy_1 = l_app(af_aw_1, Cy_rect, Af1/aw, 4);
		  n_1 = 1.5;
		  *cx_ = Round(cx_1, 3);
		  *cy_ = Round(cy_1, 3);
		  *n_ = Round(n_1, 3);
		}
	  break;
	  case CIRCULAR_TUBE:
		*cx_ = 1.26;
		*cy_ = 1.26;
		*n_ = 1.5;

	  break;
	  case DVUTAVR:
		if (index_axis == 0) { // в плоскости стенки
		  Af1 = com_sect->dvutavr.b1*com_sect->dvutavr.h1;
		  Af2 = com_sect->dvutavr.b2*com_sect->dvutavr.h2;
		  aw = com_sect->dvutavr.b*com_sect->dvutavr.h;
		  Af_max = MAX(Af1, Af2);
		  Af_min = MIN(Af1, Af2);
		  af = Af_min/Af_max;

		  cx_1 = l_app(af_aw_1, Cx_1, Af_max/aw, 4);
		  cy_1 = 1.47;
		  n_1 = 1.5;
		  if (Af_max > Af_min) {
			cx_2 = l_app(af_aw_2, Cx_2, Af_max/aw, 3);
			cy_2 = 1.47;
			n_2 = 2.0;
			*cx_ = Round(cx_1 + (1.0 - af)*(cx_2 - cx_1)/0.5, 3);
			*cy_ = Round(cy_1 + (1.0 - af)*(cy_2 - cy_1)/0.5, 3);
			*n_ = Round(n_1 + (1.0 - af)*(n_2 - n_1)/0.5, 3);
		  }
		  else {
			*cx_ = Round(cx_1, 3);
			*cy_ = Round(cy_1, 3);
			*n_ = Round(n_1, 3);
		  }
		}
		else { // в плоскости поясов
		  double cy[4] = {CY_dvutavr1};
		  aw = com_sect->dvutavr.b1*com_sect->dvutavr.h1 + com_sect->dvutavr.b2*com_sect->dvutavr.h2;
		  af = com_sect->dvutavr.b*com_sect->dvutavr.h;
		  *cx_ = 1.47;
		  *cy_ = l_app(af_aw_1, cy, af/aw, 4);;
		  *n_ = 3.0;
		}
	  break;
	  default:
		*cx_ = 1.5;
		*cy_ = 1.5;
		*n_ = 3.0;
	  }
}
//------------------------------------------------------------------------------
// Вычисление коэффициентов fi_y (7.1.3 СП), c  (п. 9.2.5 СП)
void Get_c_fi_y_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, FI_Y_PARAM *Fi_param, double fi_b, int index_comp) {
	 COMMON_SECT com_sect_comp;
	 double Af2;

	 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.h2 = thick;
	  com_sect_comp.dvutavr.b2 = Af2/thick;
	 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;
	 }

	Get_c_fi_y(DVUTAVR, &com_sect_comp, M_, N_, len_design_YZ,  length_ext, length_int,
			  param_sect, Ry, Ey, Fi_param, fi_b);
}
//------------------------------------------------------------------------------
// Вычисление коэффициентов fi_y (7.1.3 СП), c  (п. 9.2.5 СП)
void Get_c_fi_y(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, FI_Y_PARAM *Fi_param, double fi_b) {

	   int i, j, k, m, rc;
	   double lambda_uw0, lambda_uw1, lambda_uw, m_exc;
	   double tw, hw, bf, tf, hb, lambda_c, lambda_c_r, lambda_y, lambda_y_r;
	   double fi_c, fi_1, alfa_176, alfa_x, psi_, beta_, c_5, c_10, nu_;
	   double ro_, mu_, Jt_, delta_, c_max;
	   double Force[7] = {0, 0, 0, 0, 0, 0, 0};
	   double sigma_comm[4];
	   double Wc, excentr;
	   double len_ef;
	   double lambda_y_ef, lambda_y_r_ef;

	   if (type_sect != DVUTAVR)  {
		 Fi_param->fi_b = 1;
		 Fi_param->fi_y = 1;
		 Fi_param->c_ = 1;
		 return;
	   }

	   lambda_c = Round(3.14*sqrt(Ey/Ry), 3);  // стандартная гибкость  (табл. 10, обозначения)
	   lambda_c_r = 3.14;

	   fi_c = Fun_fi_713(type_sect, PI_, 1);  // значение fi_y при lambda_y=lambda_c (табл. 10, обозначения, ф. 8,9,10)

	   c_5 = 0;
	   c_10 = 0;

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

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

	   sigma_beam_comm(param_sect, Force, sigma_comm);
	   if (sigma_comm[0] < sigma_comm[2]) {
		 Wc = param_sect->Izz/param_sect->Zct;
		 if (length_ext == 0)
		   len_ef = len_design_YZ;
		 else
		   len_ef = length_ext;
	   }
	   else {
		 Wc = param_sect->Izz/(param_sect->Hsect - param_sect->Zct);
		 if (length_int == 0)
		   len_ef = len_design_YZ;
		 else
		   len_ef = length_int;
	   }
	   N_ = MAX((ABS(N_)), 1e-2);
	   if (N_ < 85000) {
		 excentr = ABS(M_/N_);
	   }
	   excentr = ABS(M_/N_);
	   m_exc = excentr*param_sect->Ar/Wc;  // ф. стр. 15 СНиП
	   m_exc = Round(m_exc, 3);

	   lambda_y = len_design_YZ/sqrt(param_sect->Iyy/param_sect->Ar);     // гибкость элемента l_efy/i
	   lambda_y_r = Round(lambda_y*sqrt(Ry/Ey), 3);  // относительная гибкость из плоскости рамы

	   lambda_y_ef = len_ef/sqrt(param_sect->Iyy/param_sect->Ar);     // гибкость элемента по длине len_ef/i
	   lambda_y_r_ef = Round(lambda_y_ef*sqrt(Ry/Ey), 3);  // относительная гибкость из плоскости рамы

		// 1. Вычисление fi_y
	   Fi_param->fi_y = Round(Fun_fi_713(type_sect, lambda_y_r, 1), 3);
	   // 2. fi_b из параметров
	   Fi_param->fi_b = fi_b;
	   // 3. Вычисление beta , табл. 10
	   //beta_ = Round(Fun_beta(lambda_y_ef, lambda_c, Fi_param->fi_y, fi_c), 2);
	   beta_ = Round(Fun_beta(lambda_y_r, lambda_c_r, Fi_param->fi_y, fi_c), 3);

		// 4. Вычисление c (ф. 57,58,59 СНиП или 9.2.4 СП 2011
	   nu_ = 1;
		//m_exc = 0.8;  Отладка
	   if (m_exc <= 5.0) {
		   alfa_x = Round(Fun_alfa(m_exc), 3);
		   //nu_ = 1 - (lambda_y_r_ef/14.0)*(2.12 - bf/param_sect->Hsect);  //ф. 97,98,99 СП
		   //Fi_param->c_ = Round(beta_/(1.0+ alfa_x*nu_*m_exc), 2);
		   Fi_param->c_ = Round(beta_/(1.0+ alfa_x*m_exc), 3);
	   }
	   else {
		   if (m_exc >= 10.0) {
			  if (Fi_param->fi_b == 0) {
				 fi_b = 0;
				 return;
			  }
			  Fi_param->c_ = 1.0/(1.0 + m_exc*Fi_param->fi_y/Fi_param->fi_b);
		   }
		   else {
			  alfa_x = Fun_alfa(5.0);
			  c_5 = beta_/(1.0+ alfa_x*5.0);
			  c_10 = 1.0/(1.0 + 10.0*Fi_param->fi_y/Fi_param->fi_b);
			  Fi_param->c_ = c_5*(2.0 - 0.2*m_exc) + c_10*(0.2*m_exc -1.0);
		   }
	   }

	   ro_=0;
		 // 4i. Вычисление с_max (ф. 60 СНиП или ф. И.1 СП)
	   ro_ = (param_sect->Izz + param_sect->Iyy)/(param_sect->Ar*hb*hb);
		 // param_sect.Itor и Jt_ отличаются множителем: 1/3 и 0.433 или 0.43
	   Jt_ = (hw*tw*tw*tw + 2*bf*tf*tf*tf)*0.43;
	   mu_ = 2.0 + 0.156*Jt_*lambda_y_ef*lambda_y_ef/(param_sect->Ar*hb*hb);
	   delta_ = 4.0*ro_/mu_;
		 // поскольку сечение симметрично, по п. И.1 СП: alfa = 0, beta = 0, B = 1.0;
		 // следовательно, ф. 60 СНиП и ф. И.1 СП совпадают
	   c_max = 2.0/(1+delta_+sqrt((1-delta_)*(1-delta_)+16*(excentr/hb)*(excentr/hb)/mu_));
	   if (lambda_y_r_ef >3.14) {   // условие п. 5.31
			Fi_param->c_ = MIN(Fi_param->c_, c_max);
	   }
	   Fi_param->c_ = Round(Fi_param->c_, 3);
	   Fi_param->c_ = MAX(Fi_param->c_, 0.001);
	   return;
}
//----------------------------------------------------------------------------
// Вычисление коэффициента fi_e для проверки устойчивости в плоскости действия момента при внецентренном сжатии
void Get_fi_e_comp(int type_sect, COMMON_SECT *com_sect, double width, double thick, double M_, double N_, int index_Y, double len_design_XZ,
				   GEOM_PARAM_SECT *param_sect, double Ry, double Ey, FI_E_PARAM *Fi_e_param, int index_comp) {
	 COMMON_SECT com_sect_comp;
	 double Af2;

	 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.h2 = thick;
	  com_sect_comp.dvutavr.b2 = Af2/thick;
	 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;
	 }

	Get_fi_e(DVUTAVR, &com_sect_comp, M_, N_,index_Y, len_design_XZ,
			  param_sect, Ry, Ey, Fi_e_param);
}
//----------------------------------------------------------------------------
// Вычисление коэффициента fi_e для проверки устойчивости в плоскости действия момента при внецентренном сжатии
void Get_fi_e(int type_sect, COMMON_SECT *com_sect, double M_, double N_, int index_Y, double len_design_XZ,
			  GEOM_PARAM_SECT *param_sect, double Ry, double Ey, FI_E_PARAM *Fi_e_param) {
	 //index_Y показывает, в какой плоскости задан момент: 0 - в плоскости стенки, 1 - в плоскости поясов
	 double m_exc, eta_, lambda_x, lambda_x_r;
	 double eta_work[4];
	 double Force[6] = {0, 0, 0, 0, 0, 0};
	 double sigma_comm[7];
	 double Wc, excentr;

	 Force[0] = N_;
	 if (index_Y==0)
	   Force[2] = M_;
	 else
	   Force[4] = M_;

	 sigma_beam_comm(param_sect, Force, sigma_comm);
	 if (index_Y == 0) {
	   if (sigma_comm[0] < sigma_comm[2])
		 Wc = param_sect->Izz/param_sect->Zct;
	   else
		 Wc = param_sect->Izz/(param_sect->Hsect - param_sect->Zct);
	 }
	 else {
		 Wc = 2*param_sect->Iyy/param_sect->Wsect;
	 }
	 N_ = MAX((ABS(N_)), 1e-2);
	 excentr = ABS(M_/N_);
	 m_exc = excentr*param_sect->Ar/Wc;  // ф. стр. 15
	 m_exc = Round(m_exc, 3);

	 if (param_sect->Izz==0 || param_sect->Iyy==0 || param_sect->Ar==0) {
		#ifdef EMBRAC
		int i = Application->MessageBox(L"Программная ошибка. Izz = 0 или Ar = 0", L" ", MB_ICONERROR);
		#else
		int i = Application->MessageBox("Программная ошибка. Izz = 0 или Ar = 0", " ", MB_ICONERROR);
		#endif
		lambda_x = len_design_XZ/100;
	 }
	 else {
	  if (index_Y == 0)
		 lambda_x = len_design_XZ/sqrt(param_sect->Izz/param_sect->Ar);     // гибкость элемента l_efx/i
	  else
		 lambda_x = len_design_XZ/sqrt(param_sect->Iyy/param_sect->Ar);     // гибкость элемента l_efx/i
	 }
	 lambda_x_r = Round(lambda_x*sqrt(Ry/Ey), 3);                       // относительная гибкость в плоскости рамы
	 // 5а. Вычисление eta (Здесь type_sect в случае составного сечения должно быть DVUTAVR, а com_sect должно учитывать наличие плиты)
	 eta_ = Fun_eta_sect(type_sect, com_sect, lambda_x_r, m_exc, eta_work, index_Y);   // табл. 73  СНиП или или Д.2 СП
	 Fi_e_param->m_ef = eta_*m_exc;
	                          // ф. 52
	 double fi_ = Fun_fi_713(type_sect, lambda_x_r, index_Y);

	 Fi_e_param->fi_e = Fun_fi_e(lambda_x_r, Fi_e_param->m_ef, fi_);          // табл. 74 СНиП или Д.3 СП

	 char stroka[120];
	 //sprintf(stroka,"N_=%6.1f кН, M_==%6.1f Кн*м, lamX=%5.3f, lamXr=%5.3f, m_exc=%5.3f, m_ef=%5.3f, eta=%5.3f, fi_e=%5.3f",
	 //				  N_*1e-3, M_*1e-6, lambda_x, lambda_x_r, m_exc, eta_, Fi_e_param->m_ef, Fi_e_param->fi_e);
	 //InfMessage("Коэфф: "+(AnsiString)stroka);

}
//-------------------------------------------------
//  Линеиная интерполяция
double lin_interpol(double x1, double x2, double y1, double y2, double x) {
   if (x1 == x2)
     return y1;
   return (y1 + (x - x1)*(y2 - y1)/(x2 - x1));
}
//-------------------------------------------------------------------------------------
// Проверка местной устойчивости стенок балок по СНиП  (п.7.2*, 7.3, 7.4*)
//                                               CП  (п.8.5.2, 8.5.3, 8.5.4)
double CheckStab_norm(SECT_DVUTAVR *sect1, SECT_DVUTAVR *sect2, double sigma_loc, double l_ef,
					  double M1, double M2, double Q_, double N_,
					  double Ry, double Rs, double Ey, double a_otsek, int index_loc, double gamma_c,
					  CHECK_STAB_PARAM *Check_stab_param ) {
// index_loc = 0 - локальная нагрузка приложена к внешнему поясу, -1 - к внутреннему

    int i, rc;
    int k, j, j1, j2, i1, i2, i0;
    double h;
    double sigma_ext, sigma_int;
	double sigma_ext1, sigma_int1;
    double sigma_ext2, sigma_int2;
    double sigma_max1, sigma_11;
    double sigma_max2, sigma_12;
    double sigma_max, tau_max, sigma_1;
    double sigma_cr, tau_cr, k_sig, k_tau;
	double krit;
    int    n_nst, j0;
    char stroka[240];
	double sigma_cr_norm_;
    double tau_cr_norm_;
    double delta_h, delta_M;
	double h0, h_ef, h_ef1, h_ef2, M_, x0, x_a, M0, N0, Q0, b_otsek, del_a;
    float  delta, c_cr;
    double a_max1, a_max2, a_max;
    double h1, h2, tw, bf, tf, hef0;
    //double lambda_w1, lambda_w2, lambda_w, lambda_ef;
    double MM;
    GEOM_PARAM_SECT param_sect;
	SECT_DVUTAVR sect0;
	double sigma_ct; // напряжение в ц.т.


    h1 = sect1->h;
	h2 = sect2->h;
	tw = sect1->b;

	Get_param_sect(sect1, &param_sect);

	sigma_ct = N_/param_sect.Ar;
	sigma_ext1 = -M1*(param_sect.Zct - sect1->h1)/param_sect.Izz + sigma_ct ;
	sigma_int1 = M1*(param_sect.Hsect - param_sect.Zct - sect2->h2)/param_sect.Izz + sigma_ct;

	// минимальное напряжение в сечении 1
	if (sigma_ext1 < sigma_int1) {
		 sigma_max1 = -sigma_ext1;  // напряжение сжатия принимается со знаком +, растяжения со знаком -
		 sigma_11 = -sigma_int1;
	}
	else {
		 sigma_max1 = -sigma_int1;  // напряжение сжатия принимается со знаком +, растяжения со знаком -
		 sigma_11 = -sigma_ext1;
	}
	if (sigma_max1 !=0.0)
	  a_max1 = Round((sigma_max1 - sigma_11)/sigma_max1, 4);
	else
	  a_max1 = 1.0;
	// расчетная высота стенки
	if (a_max1 >=1.4 && a_max1 <= 4.0)
	   h_ef1 = 2*h1/a_max1;
	else
	   h_ef1 = h1;

	//lambda_w1 = (h_ef1/tw)*sqrt(Ry/Ey);

	Get_param_sect(sect1, &param_sect);

	sigma_ct = N_/param_sect.Ar;
	sigma_ext2 = -M2*(param_sect.Zct - sect2->h1)/param_sect.Izz + sigma_ct;
	sigma_int2 = M2*(param_sect.Hsect - param_sect.Zct - sect2->h2)/param_sect.Izz + sigma_ct;

	// минимальное напряжение в сечении 1
	if (sigma_ext2 < sigma_int2) {
		 sigma_max2 = -sigma_ext2;  // напряжение сжатия принимается со знаком +, растяжения со знаком -
		 sigma_12 = -sigma_int2;
	}
	else {
		 sigma_max2 = -sigma_int2;  // напряжение сжатия принимается со знаком +, растяжения со знаком -
		 sigma_12 = -sigma_ext2;
	}
	if (sigma_max2!=0.0)
	  a_max2 = Round((sigma_max2 - sigma_12)/sigma_max2, 4);
	else
	  a_max2 = 1;
	// расчетная высота стенки
	if (a_max2 >=1.4 && a_max2 <= 4.0)
	   h_ef2 = 2*h1/a_max2;
	else
	   h_ef2 = h2;

	//lambda_w2 = (h_ef2/tw)*sqrt(Ry/Ey);

	hef0 = (h_ef1 + h_ef2)/2;

	// определение расчетного сечения
	MM = MAX((ABS(M1)), (ABS(M2)));
	if (MM>0)
	  M_ = ABS((M1 - M2)/MM);
	else
	  M_ = 0;
	b_otsek = 0;
	x0 = 0;
	if (M_ <= 1.0 ) {  // момент не меняет знак
	   if (a_otsek <= hef0) {
		  x_a = a_otsek/2;
	   }
	   else {
		  if ((ABS(M1)) > (ABS(M2)))
			   x_a = hef0/2;
			else
			   x_a = a_otsek - hef0/2;
	   }
	 }
	 else {              // момент меняет знак
		 if ((ABS(M1)) > (ABS(M2))) {
			b_otsek = a_otsek*(-M1)/(-M1+M2);
			if (b_otsek <= hef0) {
			   x_a = b_otsek/2;
			}
			else {
			   x_a = hef0/2;
			}
		 }
		 else  {
			x0 = a_otsek*(-M1)/(-M1+M2);
			b_otsek = a_otsek - x0;
			if (b_otsek <= hef0) {
			   x_a = x0 + b_otsek/2;
			}
			else {
			   x_a = a_otsek - hef0/2;
			}
		 }
	   }

	   delta_h = (h2-h1)/a_otsek;
	   delta_M = (M2-M1)/a_otsek;
	   //
	   //sprintf(stroka,"delta_M=%g, x_a=%g",  delta_M,  x_a);
	   //InfMessage(stroka);
	   //
	   M0 = M1 + delta_M*x_a;
	   h0 = h1 + delta_h*x_a;

	   memcpy(&sect0, sect1, sizeof(SECT_DVUTAVR));
	   sect0.h = h0;
	   Get_param_sect(&sect0, &param_sect);

	   sigma_ct = N_/param_sect.Ar;
	   sigma_ext = -M0*(param_sect.Zct - sect0.h1)/param_sect.Izz + sigma_ct;
	   sigma_int = M0*(param_sect.Hsect - param_sect.Zct - sect0.h2)/param_sect.Izz + sigma_ct;

       if (sigma_ext < sigma_int) {
             sigma_max = -sigma_ext;  // напряжение сжатия принимается со знаком +, растяжения со знаком -
             sigma_1 = -sigma_int;
             bf = sect1->b1;
             tf = sect1->h1;
       }
       else {
             sigma_max = -sigma_int;  // напряжение сжатия принимается со знаком +, растяжения со знаком -
             sigma_1 = -sigma_ext;
			 bf = sect1->b2;
             tf = sect1->h2;
       }

       if (sigma_max <= 0) {
         return 0.0;
       }
	   a_max = Round((sigma_max - sigma_1)/sigma_max, 4);
       //a_max = MAX(1, a_max);
       // a_max=2.0;
       // приведенная высота стенки, необходима для функций sigma_cr_norm(), tau_cr_norm_ = tau_cr_norm(a_otsek[i]);
       if (a_max >=1.4 && a_max <= 4.0)
         h_ef = 2*h0/a_max;
	   else
         h_ef = h0;
       //lambda_w = (h_ef/tw)*sqrt(Ry/Ey);
	   //lambda_ef = ((MIN(h_ef, a_otsek)))*sqrt(Ry/Ey)/tw;

       tau_max = (ABS(Q_))/(tw*h0);

       sigma_cr = Sigma_cr_norm(Ry, Ey, h_ef, bf, tf, tw, &delta, &c_cr);
       tau_cr = Tau_cr_norm(Ry, Rs, Ey, h_ef, tw, a_otsek);
	   // Обработка локальных напряжений
	   double beta = 0.8;
	   double sigma_loc_cr = Sigma_loc_cr(Ry, Ey, h_ef, l_ef, beta, bf, tf, tw, a_otsek);

	   k_sig = sigma_max/sigma_cr;
	   k_tau = tau_max/tau_cr;
	   double k_sig_loc = sigma_loc/sigma_loc_cr;

       //sprintf(stroka,"sigma_cr=%g, tau_cr=%g, a_max=%g",  sigma_cr, tau_cr, a_max);
	   //InfMessage(stroka);
	   //sprintf(stroka,"sigma=%g, tau=%g", sigma_max, tau_max);
	   //InfMessage(stroka);
	   if (sigma_ext < sigma_int && index_loc == 0 ||
		   sigma_ext > sigma_int && index_loc == 1) {
		   krit = sqrt((k_sig + k_sig_loc)*(k_sig + k_sig_loc) + k_tau*k_tau)/gamma_c;
	   }
	   else {
		   krit = MAX(sqrt((k_sig)*(k_sig) + k_tau*k_tau)/gamma_c,
					  sqrt((k_sig_loc)*(k_sig_loc) + k_tau*k_tau)/gamma_c);
	   }
	   Check_stab_param->M1 = M1;
	   Check_stab_param->M2 = M2;
	   Check_stab_param->Q_ = Q_;
	   Check_stab_param->N_ = N_;
	   Check_stab_param->sigma_max = sigma_max;
	   Check_stab_param->sigma_cr = sigma_cr;
	   Check_stab_param->tau_max = tau_max;
	   Check_stab_param->tau_cr = tau_cr;
	   Check_stab_param->krit = krit;
	   Check_stab_param->sigma_loc = sigma_loc;
	   Check_stab_param->sigma_loc_cr = sigma_loc_cr;
	   return krit;
}
//------------------------------------------------------
// функции определения критических напряжений потери устойчивости
double Sigma_cr_norm(double Ry, double Ey, double hef, double bf, double tf, double tw, float *delta_, float* c_cr_) {
    double sig_, lambda_w_max_, c_cr, kR;
    kR = sqrt(Ry/Ey);
    lambda_w_max_ = hef*kR/tw;
    c_cr = cr(hef, bf, tf, tw, delta_);
    *c_cr_ = c_cr;
    sig_ = c_cr*Ry/(lambda_w_max_*lambda_w_max_);
    return sig_;
}
//----------------------------------------------------------------
double Sigma_loc_cr(double Ry, double Ey, double hef, double lef, double beta, double bf, double tf, double tw, double a_otsek) {
	double sig_, lambda_w_max_, c1, c2, kR;
    double ro_, delta_;
	double Koeff_I_J;
    kR = sqrt(Ry/Ey);

    ro_ = 1.04*lef/hef;
	delta_ = beta*(bf/hef)*pow(tf/tw, 3);    // ф. 84
	c1 = App_matrix_fun(Tabl_14[0], 7, 10, Column_14, Line_14, a_otsek/hef, ro_);
    //App_matrix(Tabl_15, Line_15, Column_15, 6, 8, a_otsek/hef, delta_, 1.6, 30.0);
    c2 = App_matrix_fun(Tabl_15[0], 6, 8, Column_15, Line_15, a_otsek/hef, delta_);

    lambda_w_max_ = hef*kR/tw;
    sig_ = c1*c2*Ry/(lambda_w_max_*lambda_w_max_);
    return sig_;
}
//-------------------------------------------------------
double Tau_cr_norm(double Ry, double Rs, double Ey, double hef, double tw, double a_otsek) {
    double tau_, mu, lambda_ef_max_, kR;
    kR = sqrt(Ry/Ey);
    lambda_ef_max_ = (MIN(hef, a_otsek))*kR/tw;
    if (lambda_ef_max_==0) return 1e6;
    mu = a_otsek/hef;
    if (mu<1) mu =1/mu;
    tau_ = 10.3*(1+0.76/(mu*mu))*(Rs/(lambda_ef_max_*lambda_ef_max_));
    return tau_;
}

//-------------------------------------------------------
double cr(double hef, double bf, double tf, double tw, float *delta_) {
   double delta, c_;

   delta = 0.8*(bf/hef)*pow(tf/tw, 3);
   if (delta<0.8) return 30.0;
   if (delta>30) return 35.5;
   *delta_ = delta;
   c_ = l_app(delta_cr, c_cr, delta, 7);
   return c_;
}
//--------------------------------------------------------------
// Сравнение сечений
int sect_cmp(DATA_SECT_SCAD *sect1, DATA_SECT_SCAD *sect2, int index_cmp) {
  if (sect1->type_sect != sect2->type_sect)
    return 1;
  if (sect1->com_sect.dvutavr.n_group >=0 ) {   // стандартный профиль
    if (sect1->com_sect.dvutavr.n_group != sect2->com_sect.dvutavr.n_group)
      return 2;
    if (sect1->com_sect.dvutavr.n_profil != sect2->com_sect.dvutavr.n_profil)
      return 2;
  }
  else {                       // сварной профиль
    if (sect1->com_sect.dvutavr.b != sect2->com_sect.dvutavr.b)
      return 3;
    if (index_cmp==0) {
      if (sect1->com_sect.dvutavr.h != sect2->com_sect.dvutavr.h)
        return 4;
    }
    if (sect1->com_sect.dvutavr.b1 != sect2->com_sect.dvutavr.b1)
      return 5;
    if (sect1->com_sect.dvutavr.b2 != sect2->com_sect.dvutavr.b2)
      return 5;
    if (sect1->com_sect.dvutavr.h1 != sect2->com_sect.dvutavr.h1)
      return 5;
    if (sect1->com_sect.dvutavr.h2 != sect2->com_sect.dvutavr.h2)
      return 5;
    if (sect1->com_sect.dvutavr.n_group != sect2->com_sect.dvutavr.n_group)
      return 5;
    if (sect1->com_sect.dvutavr.n_profil != sect2->com_sect.dvutavr.n_profil)
      return 5;
  }
  return 0;
}
//--------------------------------------------------------------
// Сравнение сечений двутавра
int sect_dvut_cmp(SECT_DVUTAVR *sect1, SECT_DVUTAVR *sect2, int index_cmp) {
  if (sect1->n_group >=0 ) {   // стандартный профиль
    if (sect1->n_group != sect2->n_group)
      return 2;
    if (sect1->n_profil != sect2->n_profil)
      return 2;
  }
  else {                       // сварной профиль
    if (sect1->b != sect2->b)
      return 3;
    if (index_cmp==0) {
      if (sect1->h != sect2->h)
        return 4;
    }
    if (sect1->b1 != sect2->b1)
      return 5;
    if (sect1->b2 != sect2->b2)
      return 5;
    if (sect1->h1 != sect2->h1)
      return 5;
    if (sect1->h2 != sect2->h2)
      return 5;
    if (sect1->n_group != sect2->n_group)
      return 5;
    if (sect1->n_profil != sect2->n_profil)
      return 5;
  }
  return 0;
}
//--------------------------------------------------------------
// Сравнение данных о сечениях
int data_sect_cmp(DATA_SECT *sect1, DATA_SECT *sect2, int index_cmp) {
  if (sect1->local_sect != sect2->local_sect)
    return 1;
  if (sect1->pos_sect != sect2->pos_sect)
    return 2;

  if (sect1->section.n_group >=0 ) {   // стандартный профиль
    if (sect1->section.n_group != sect2->section.n_group)
      return 3;
    if (sect1->section.n_profil != sect2->section.n_profil)
      return 3;
  }
  else {                       // сварной профиль
    if (sect1->section.b != sect2->section.b)
      return 4;
    if (index_cmp==0) {
      if (sect1->section.h != sect2->section.h)
        return 5;
    }
    if (sect1->section.b1 != sect2->section.b1)
      return 6;
    if (sect1->section.b2 != sect2->section.b2)
      return 6;
    if (sect1->section.h1 != sect2->section.h1)
      return 6;
    if (sect1->section.h2 != sect2->section.h2)
      return 6;
    if (sect1->section.n_group != sect2->section.n_group)
      return 6;
    if (sect1->section.n_profil != sect2->section.n_profil)
      return 6;
  }
  return 0;
}
//--------------------------------------------------------------
// Сравнение данных о ребрах
int data_edge_cmp(EDGE_PARAM *sect1, EDGE_PARAM *sect2, int index_cmp) {
  if (sect1->local_sect != sect2->local_sect ||
      sect1->pos_common != sect2->pos_common ||
      sect1->pos_single != sect2->pos_single ||
      sect1->pos_sect != sect2->pos_sect ||
      sect1->flag_double != sect2->flag_double ||
      sect1->flag_all != sect2->flag_all ||
      sect1->b_edge!= sect2->b_edge ||
      sect1->t_edge!= sect2->t_edge ||
      sect1->h_edge!= sect2->h_edge ||
      sect1->pos_edge_norm!= sect2->pos_edge_norm ||
      sect1->t_pop!= sect2->t_pop)
    return 1;

    if (index_cmp==0) {
      if (sect1->h_stenki != sect2->h_stenki)
        return 2;
    }
  return 0;
}
//--------------------------------------------------------------
// Сравнение данных о продольных ребрах
int data_long_edge_cmp(LONG_EDGE_PARAM *sect1, LONG_EDGE_PARAM *sect2) {
  if (sect1->local_sect != sect2->local_sect ||
      sect1->pos_polka != sect2->pos_polka ||
      sect1->index_polka != sect2->index_polka ||
      sect1->index_single != sect2->index_single ||
      sect1->index_pos_left != sect2->index_pos_left ||
      sect1->index_pos_right != sect2->index_pos_right ||
      sect1->pos_sect1 != sect2->pos_sect1 ||
      sect1->pos_sect2 != sect2->pos_sect2 ||
      sect1->flag_double != sect2->flag_double ||
      sect1->flag_all != sect2->flag_all ||
      sect1->b_edge!= sect2->b_edge ||
      sect1->t_edge!= sect2->t_edge)
    return 1;

  return 0;
}
//-------------------------------------------------------
//  Функция аппроксимации по матрице
double App_matrix_fun(double *Tabl_0, int nrow, int ncol, double *Column, double *Row, double x_, double y_) {
// nrow - число сткок матрицы
// ncol - число столбцов матрицы
// Tabl - матрица ncol*nrow
// Column - значения ординат (длина nrow)
// Row - значения абсцисс (длина ncol)
// x_ - абсцисса функции
// y_ - ордината функции
   double value, di, dj, g1, g2;
   int i1, j1;
   double **Tabl;
   int rc;

   ALLOCMAT(Tabl, nrow, ncol);              //Tabl[0]
   //Blpmat(nrow, ncol, 0, Tabl_0, Tabl);

			for (i1=0; i1<nrow; i1++)
			  if (Column[i1] > y_) break;
			i1--;
			for (j1=0; j1<ncol; j1++)
			  if (Row[j1]>x_) break;
			j1--;
			i1 = ((nrow-1)<(i1)) ? (nrow-1):(i1);
			j1 = ((ncol-1)<(j1)) ? (ncol-1):(j1);
			if (i1 == -1) {
			  di = 0.0;
			  i1=0;
			}
			else {
			  if (i1==nrow-1)
				di = 0;
			  else
				di = (y_- Column[i1])/(Column[i1+1] - Column[i1]);
			}
			if (j1==-1) {
			  j1 = 0;
			  dj = 0;
			}
			else {
			  if (j1==ncol-1)
				dj = 0;
			  else
				dj = (x_ - Row[j1])/(Row[j1+1] - Row[j1]);
			}
			int i12 = i1+1;
			int i22 = i1+1;
			int j21 = j1+1;
			int j22 = j1+1;
			if (i1==nrow-1) {
			  i12 = i1;
			  i22 = i1;
			}
			if (j1==ncol-1) {
			  j21 = j1;
			  j22 = j1;
			}


			//g1 = Tabl[i1][j1]*(1 - di) + Tabl[i12][j1]*di;
			//g2 = Tabl[i1][j21]*(1 - di) + Tabl[i22][j22]*di;
			g1 = Tabl_0[i1*ncol + j1]*(1 - di) + Tabl_0[i12*ncol + j1]*di;
			g2 = Tabl_0[i1*ncol + j21]*(1 - di) + Tabl_0[i22*ncol + j22]*di;

			value = g1*(1 - dj) + g2*dj;
			FREEMAT(Tabl);
			return value;
}
#include "ElementCheckFunction%Stab_teta.cpp"

#pragma package(smart_init)









