Unit Refinar;
Interface

Uses Crt, AccPampa, gnuplot, StrUtils, Math,SysUtils, Classes, FileUtil, Forms, Controls, Graphics, Dialogs, ValEdit,
ExtCtrls, Menus, ComCtrls, StdCtrls;


{Procedimientos}

Procedure Leer; {lee el archivo .dat y guarda toda la informacion}
Procedure Leer_esp; {lee el archivo .spe y guarda toda la informacion}
Procedure Escribir; {escribe el nuevo .dat }
Procedure EscribirSal; {escribe el nuevo .dat }
Procedure Graficar; {grafica el espectro experimental y el de pampa}
Procedure Archivos; {lee los archivos .asc (contienen informacion que usa el programa) y guarda la informacion}
Procedure Refina; {Refina los parametros}
Procedure Picosuma; {Se calcula la intensidad de los picos sumas y la contribucion que hace cada pico a estos}
Procedure ContNPR(); {Cuenta los parametros que se refinan}
Procedure TablaEfi(); {Escribe un archivo con la eficiencia y energía}
Procedure Errores;
Procedure TablaSP(); {Llena el arreeglo de SpAr, leyendo los archivos del stopping del .dat}


{Funciones}

Function Perfil_Gaus(E:Double;Var En:Double;sigma2:Double) :Double;  {Forma gaussiana del pico}
Function Perfil_RAS(E: Double;Var En: Double;sigma2: Double) :Double;
Function Intensidad(U:Integer;J:Integer;E:Double;S:Double) :Double;  {Intensidad del pico}
Function Intensidad_film(U:Integer;J:Integer;E:Double;S:Double) :Double;  {Intensidad del pico del film}
Function Intensidad_mg(U:Integer;J:Integer;E:Double) :Double;  {Intensidad del pico par muestras gruesas}
Function Intensidad_mg2(U:Integer;J:Integer;E:Double) :Double;  {Intensidad del pico par muestras gruesas}
Function Mu(Z:Integer;E:Double) :Double;  {coeficiente de atenuacion}
Function MuM(E:Double):Double; {Mu de la muestra (se usa para muestras gruesas)}
Function MuM2(E:Double):Double; {promedio de los Mu multiplicados por la densidad}
Function Efi(Ei:Double) :Double; {Calculo de la eficiencia del detector}
Function Fondo(E:Double;A:Double;B:Double;C:Double) :Double; {Fondo bremsstrahlung}
Function FondoProm(I:Integer;d:Double):Double;   {Fondo para muestras gruesas(promedio del fondo en 3sigma para suavisar la atenuacion)}
Function FondAt(En:Double;depth:Double) :Double;    {atenuacion del fondo para muestras gruesas}
Function Func(P :DoubleArrayNP; NPR :Double) :Double;  {Calcula la intensidad de cada canal y entrega el chi^2}
Function Amoeba(PPP:pMN; NPR:Integer) :DoubleArrayNP;  {Entrega un conjunto de parametros minimizado}
Function SigmaK_ion(Z:Integer) :Double; {Seccion eficaz de ionizacion: polinomio de grado 6 ajustado a los valores de Lappicki}
Function SigmaK(Z:Integer;Ei:Double) :Double; {datos de PIXE2010: Proton/alpha Ionization (k, l, m shell), Tabulated Cross Section Library.}
Function SigmaL1(Z:Integer;Ei:Double) :Double; {datos de PIXE2010: Proton/alpha Ionization (k, l, m shell), Tabulated Cross Section Library.}
Function SigmaL2(Z:Integer;Ei:Double) :Double; {datos de PIXE2010: Proton/alpha Ionization (k, l, m shell), Tabulated Cross Section Library.}
Function SigmaL3(Z:Integer;Ei:Double) :Double; {datos de PIXE2010: Proton/alpha Ionization (k, l, m shell), Tabulated Cross Section Library.}
Function SigmaM(Z:Integer;Ei:Double) :Double; {datos de PIXE2010: Proton/alpha Ionization (k, l, m shell), Tabulated Cross Section Library.}
Function SPM(Ei:Double) :Double; {stopping power medio con datos de stopping power sacados de SRIM.}
Function Resc(U,J:Integer):Double; {Pico de escape}
Function Trans(k:Integer;Ener:Double;u,j:Integer) :Double; {Coeficiente de transmision para muestras gruesas}
Function Trans2(k:Integer;Ener:Double;u,j:Integer) :Double; {Coeficiente de transmision para muestras gruesas}
Function RhoZmax() :Double; {calculo de la profundidad maxima en una muestra gruesa}


Implementation
uses Unit1;
Var
   P_final:doublearraynp;

Procedure Leer;
    Var  
    Dat,Esp        :Text; {Para leer el archivo .dat y .spe} 
    ch             :Char;
    i,j,l          :Integer; {para los for} 
    Data           :String[6];
    Dato           :string[14];
    DiEsp2         :String;
Begin
  
    Assign(Dat,DiDat);
    Reset(Dat);
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,ArEsp);
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,comentario);
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,DiEsp);
    Repeat Read(Dat,Ch) Until Ch=':';
    Read(Dat,Einc);
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,Detector);
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,Toff); {angulo de take off}
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,Nel);  {numero total de elementos}
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,Ndet); {numero de elementos detectados}
    Nnodet:=Nel-Ndet;
    Repeat Read(Dat,Ch) Until Ch=':';
    For i:=1 To Ndet Do Begin   {numeros atomicos de los detectados}
        Read(Dat,Zdet[i]);
    End;
    Repeat Read(Dat,Ch) Until Ch=':';
    For i:=1 To Nnodet Do Begin {numeros atomicos de los no detectados}
        Read(Dat,Znodet[i]);
    End;
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,Nzonas);
    Readln(Dat,Canali,Canalf); {canal inicial y final de la zona de interes}
    Readln(Dat); Readln(Dat);
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,P0[1],Ch,Ch,Ch,PE[1],Cod[1],Lam[1]);         {parametro 1 = ganancia}
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,P0[2],Ch,Ch,Ch,PE[2],Cod[2],Lam[2]);         {parametro 2 = cero}
    Readln(Dat); Readln(Dat);
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,P0[3],Ch,Ch,Ch,PE[3],Cod[3],Lam[3]);         {parametro 3 = espesor de Be en cm}
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,P0[4],Ch,Ch,Ch,PE[4],Cod[4],Lam[4]);         {parametro 4 = espesor de Au o Ni en cm}
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,P0[5],Ch,Ch,Ch,PE[5],Cod[5],Lam[5]);         {parametro 5 = espesor de CM(Si) en cm}
    Readln(Dat); Readln(Dat);
    Readln(Dat,P0[6],Ch,Ch,Ch,PE[6],Cod[6],Lam[6]);         {parametro 6 = escala del fondo}
    Readln(Dat,P0[7],Ch,Ch,Ch,PE[7],Cod[7],Lam[7]);         {parametro 7 = escala del fondo}
    Readln(Dat,P0[8],Ch,Ch,Ch,PE[8],Cod[8],Lam[8]);         {parametro 8 = escala del fondo}
    Readln(Dat); Readln(Dat); Readln(Dat);
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,P0[9],Ch,Ch,Ch,PE[9],Cod[9],Lam[9]);         {parametro 9 = Noise}
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,P0[10],Ch,Ch,Ch,PE[10],Cod[10],Lam[10]);         {parametro 10 = Factor de Fano}
    Readln(Dat); Readln(Dat);
    NL:=0;
    For i:=1 To NDet Do Begin	           {aca se lee el numero de lineas para cada elemento y se suman todos en NL}
        Repeat Read(Dat,Ch) Until Ch=':';
        Readln(Dat,Nlin[i]);
        NL:=NL+Nlin[i];
        For j:=1 To Nlin[i]+1 Do Readln(Dat)
    End;
    Reset(Dat);       {Se lee desde el principio el archivo}
    For i:=1 To 33 Do Readln(Dat);
    l:=NPG;
    For i:=1 To Ndet Do Begin
        For j:=1 To Nlin[i] do Begin
            Inc(l);
            Read(Dat,Perf[i,j],P0[l],Ch,Ch,Ch,PE[l],P0[NL+l],Ch,Ch,Ch,PE[NL+l],P0[2*NL+l],Ch,Ch,Ch,PE[2*NL+l],P0[3*NL+l],Ch,Ch,Ch,PE[3*NL+l]);
            {parametros 13 a 12+NL = energias , 13+NL a 12+2*NL = fracciones, 13+2*NL a 12+3*NL = tail, 13+3*NL a 12+4*NL = Beta}																		 {parametros 10+Nl al 9+2*NL = fracciones de linea}
	    Repeat Read(Dat,Borde[i,j][1]) Until Borde[i,j][1]<>' ';
	    Read(Dat,Borde[i,j][2]);
	    Repeat Read(Dat,Desc[l-NPG]) Until Desc[l-NPG]<>' ';
	    Readln(Dat,Cod[l],Cod[NL+l],Cod[2*NL+l],Cod[3*NL+l],Lam[l]);
	    Perf[i,j]:=UpCase(Perf[i,j]); {transforma a mayuscula}
            If (Perf[i,j]<>'V') And (Perf[i,j]<>'G') Then Perf[i,j]:='V';
            Lam[NL+l]:=Lam[l];
            Lam[2*NL+l]:=Lam[l];
            Lam[3*NL+l]:=Lam[l];
            E[i,j]:=P0[l];      {Guarda las energias}
	    Fr[i,j]:=P0[NL+l];  {Guarda las fracciones}
            Tail[i,j]:=P0[2*NL+l]; {Guarda las colas}
            Beta[i,j]:=P0[3*NL+l]; {Guarda los beta}
	End;
	Readln(Dat);
	Readln(Dat)
    End;
    NG:=-1;
    Repeat
    Readln(Dat,Data);
    Inc(NG)
    Until Data='';
    Reset(Dat);
    Repeat Readln(Dat,Data) Until Data='Grupos';
    For i:=1 To NG Do Begin         							{parametros 13+4*NL al 12+4*NL+NG = Factor de escala de grupo}
	Readln(Dat,Grupo[i],P0[NPG+4*NL+i],Ch,Ch,Ch,PE[NPG+4*NL+i],Cod[NPG+4*NL+i],Lam[NPG+4*NL+i]); {Guarda los grupos}
	Fesc[i]:=P0[NPG+4*NL+i];  {**Factores de escala**Deben estar en este orden: KK, L3, L2, L1, M5, M4, M3, ...}
    End;
    ReadLn(Dat); ReadLn(Dat);
    For i:=NPG+NG+4*NL+1 To NPG+NG+4*NL+Ndet Do Begin          {parametros 13+4*NL+NG al 12+4*NL+NG+Nel = concentraciones}
	Read(Dat,Zdet[i-NPG-NG-4*NL],P0[i],Ch,Ch,Ch,PE[i],Cod[i],Lam[i]);
	C[i-(NPG+NG+4*Nl)]:=P0[i];  {**Concentraciones**En el orden en que estan en el .dat}
    End;
    For i:=NPG+NG+4*NL+Ndet+1 To NPG+NG+4*NL+Nel Do Begin
	Read(Dat,Znodet[i-NPG-NG-4*NL-Ndet],P0[i],Ch,Ch,Ch,PE[i],Cod[i],Lam[i]);
	C[i-(NPG+NG+4*Nl)]:=P0[i];  {**Concentraciones**En el orden en que estan en el .dat}
    End;
    {**arreglo Z con los numeros atomicos de todos los elementos**}
    For i:=1 To Ndet Do Z[i]:=Zdet[i];
    For i:=Ndet+1 To Nel Do Z[i]:=Znodet[i-Ndet];
    Readln(Dat);
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,Titer);     {Lee la Tolerancia para las iteraciones}
    Readln(Dat);
    Readln(Dat);			{parametro 11 a 13= constante pico suma}
    Readln(Dat,P0[11],Ch,Ch,Ch,PE[11],Cod[11],Lam[11]);
    Readln(Dat,P0[12],Ch,Ch,Ch,PE[12],Cod[12],Lam[12]);
    Readln(Dat,P0[13],Ch,Ch,Ch,PE[13],Cod[13],Lam[13]);
    Readln(Dat);
    Readln(Dat);
    Readln(Dat,P0[14],Ch,Ch,Ch,PE[14],Cod[14],Lam[14]); {parametro 14 = intensidad pico de fluorescencia interna}
    Readln(Dat);
    Repeat Read(Dat,Ch) Until Ch=':';
    Readln(Dat,Espesor);                   {Lee el espesor}
    {Comenzamos a leer el film}
    Nfilm:=0;
    Nfdet:=0;
    Nfnodet:=0;
    NLF:=0;
    If Not SeekEof(Dat) Then Begin
       Readln(Dat);
       Repeat Read(Dat,Ch) Until Ch=':';
       Readln(Dat,Nfilm);      {Lee el numero de elementos en el film}
       Repeat Read(Dat,Ch) Until Ch=':';
       Readln(Dat,Nfdet);      {Lee el numero de elementos detectables en el film}
       Nfnodet:=Nfilm-Nfdet;
       Repeat Read(Dat,Ch) Until Ch=':';
       For i:=Ndet+1 To Ndet+Nfdet Do Begin {lee los Z de los elementos detectables del film}
           Read(Dat,Zdet[i]);
       End;
       Repeat Read(Dat,Ch) Until Ch=':';
       For i:=Nnodet+1 To Nnodet+Nfnodet Do Begin {lee los Z de los elementos no detectables del film}
            Read(Dat,Znodet[i]);
       End;
       Readln(Dat);Readln(Dat);Readln(Dat);
       For i:=Ndet+1 To Ndet+Nfdet Do Begin	    {aca se lee el numero de lineas para cada elemento y se suman todos en NLF}
            Repeat Read(Dat,Ch) Until Ch=':';
	    Readln(Dat,Nlin[i]);
	    NLF:=NLF+Nlin[i];
	    For j:=1 To Nlin[i]+1 Do Readln(Dat)
       End;
       Reset(Dat);       {Se lee desde el principio el archivo}
       Repeat Readln(Dat,Dato) Until Dato='Picos del film';
       Readln(Dat);
       l:=NPG;
       For i:=Ndet+1 To Ndet+Nfdet Do Begin
            For j:=1 To Nlin[i] do Begin
	      Inc(l);
	      Read(Dat,Perf[i,j],P0[4*NL+NG+Nel+l],Ch,Ch,Ch,PE[4*NL+NG+Nel+l],P0[4*NL+NG+Nel+NLF+l],Ch,Ch,Ch,PE[4*NL+NG+Nel+NLF+l],
              P0[4*NL+NG+Nel+2*NLF+l],Ch,Ch,Ch,PE[4*NL+NG+Nel+2*NLF+l],P0[4*NL+NG+Nel+3*NLF+l],Ch,Ch,Ch,PE[4*NL+NG+Nel+3*NLF+l]);
	      {parametros 11+2*NL+NG+Nel a 10+2*NL+NG+Nel+NLF = energias}
	      {parametros 11+2*NL+NG+Nel+NLF al 10+2*NL+NG+Nel+2*NLF = fracciones de linea}
              {parametros 11+2*NL+NG+Nel+2*NLF al 10+2*NL+NG+Nel+3*NLF = colas}
              {parametros 11+2*NL+NG+Nel+3*NLF al 10+2*NL+NG+Nel+4*NLF = betas}
	      Repeat Read(Dat,Borde[i,j][1]) Until Borde[i,j][1]<>' ';
	      Read(Dat,Borde[i,j][2]);
	      Repeat Read(Dat,Desc[NL+l-NPG]) Until Desc[NL+l-NPG]<>' ';
	      Readln(Dat,Cod[4*NL+NG+Nel+l],Cod[4*NL+NG+Nel+NLF+l],Cod[4*NL+NG+Nel+2*NLF+l],Cod[4*NL+NG+Nel+3*NLF+l],Lam[4*NL+NG+Nel+l]);
	      Perf[i,j]:=UpCase(Perf[i,j]); {transforma a mayuscula}
	      If (Perf[i,j]<>'V') And (Perf[i,j]<>'G') Then Perf[i,j]:='V';
	      Lam[4*NL+NG+Nel+NLF+l]:=Lam[4*NL+NG+Nel+l];
              Lam[4*NL+NG+Nel+2*NLF+l]:=Lam[4*NL+NG+Nel+l];
              Lam[4*NL+NG+Nel+3*NLF+l]:=Lam[4*NL+NG+Nel+l];
	      E[i,j]:=P0[4*NL+NG+Nel+l];      {Guarda las energias}
	      Fr[i,j]:=P0[4*NL+NG+Nel+NLF+l];  {Guarda las fracciones}
	      Tail[i,j]:=P0[4*NL+NG+Nel+2*NLF+l];  {Guarda las colas}
	      Beta[i,j]:=P0[4*NL+NG+Nel+3*NLF+l];  {Guarda los betas}
            End;
	    Readln(Dat);
	    Readln(Dat)
       End;
       l:=NPG;
       For i:=Ndet+1 To Ndet+Nfdet Do Begin          {parametros 11+2*NL+NG+Nel+2*NLF al 10+2*NL+NG+Nel+2*NLF+Nfilm = concentraciones}
            Inc(l);
            Read(Dat,Zdet[i],P0[4*NL+NG+Nel+4*NLF+l],Ch,Ch,Ch,PE[4*NL+NG+Nel+4*NLF+l],Cod[4*NL+NG+Nel+4*NLF+l],Lam[4*NL+NG+Nel+4*NLF+l]);
	    C[i-Ndet+Nel]:=P0[4*NL+NG+Nel+2*NLF+l];  {**Concentraciones**En el orden en que estan en el .dat}
       End;
       For i:=Nfdet+1 To Nfilm Do Begin
           Inc(l);
           Read(Dat,Znodet[i-Nfdet+Nnodet],P0[4*NL+NG+Nel+4*NLF+l],Ch,Ch,Ch,PE[4*NL+NG+Nel+4*NLF+l],Cod[4*NL+NG+Nel+4*NLF+l],Lam[4*NL+NG+Nel+4*NLF+l]);
	   C[i+Nel]:=P0[4*NL+NG+Nel+4*NLF+l];  {**Concentraciones**En el orden en que estan en el .dat}
       End;
       {**arreglo Z con los numeros atomicos de todos los elementos**}
       For i:=Nel+1 To Nel+Nfdet Do Z[i]:=Zdet[i-Nnodet];
       For i:=Nel+Nfdet+1 To Nel+Nfilm Do Z[i]:=Znodet[i-Ndet-Nfdet];
       Repeat Read(Dat,Ch) Until Ch=':';
       Readln(Dat,Espesorfilm);                   {Lee el espesor}
    End;
    ContNPR;
    Close(Dat);
End;

Procedure Leer_Esp;
    Var
    Esp            :Text; {Para leer el archivo .dat y .spe}
    i              :Integer; {para los for}
    DiEsp2         :String;
Begin
    Case SistOp of         {Lee el espectro esperimental}
         1: DiEsp2 := DiEsp+'/';
         2: DiEsp2 := DiEsp+'\';
    End;
    Assign(Esp,DiEsp2+ArEsp);
    Reset(Esp);
    Readln(Esp);
    Readln(Esp,i,NCanal);
    NCanal:=NCanal+1-i;
    For i:=1 to 30000 Do Ymed[i]:=0;
    For i:=1 To NCanal Do Read(Esp,Ymed[i]);
    Close(Esp);

End;  

Procedure ContNPR();
Var
    i:Integer;
Begin
    For i:=1 To NPG+4*NL+NG+Nel+4*NLF+Nfilm Do Begin
        Sumcod[i]:=Sumcod[i-1]+Cod[i];
        NPR:=Sumcod[i];{numero de parametros a refinar}
    End;
    Form1.label6.Caption:=FloatToStr(NPR);
    Form1.Refresh;
End;

Procedure Escribir;
    Var  
    Dat            :Text; {Para leer el archivo .dat y .spe} 
    i,j,l          :Integer; {para los for} 
Begin

    Assign(Dat,DiDat2);
    Rewrite(Dat);
    Writeln(Dat,'Espectro:'+ArEsp);
    Writeln(Dat,'Comentario:'+comentario);
    Writeln(Dat,'Directorio:'+DiEsp);
    Writeln(Dat,'Energía de incidencia (keV):',Einc:4:1,'      Espectrómetro:',Detector);
    Writeln(Dat,'Angulo de salida (grados): ',Toff:4:1);
    Writeln(Dat,'Número total de elementos: ',Nel:2);
    Writeln(Dat,'Número de elementos detectados: ',Ndet:2);
    Write(Dat,'Números atómicos de los detectables: ');
    For i:=1 To Ndet Do Write(Dat,Zdet[i]:3);
    Writeln(Dat);
    Write(Dat,'Números atómicos de los no detectables: ');
    For i:=1 To Nnodet Do Write(Dat,Znodet[i]:3);
    Writeln(Dat);
    Writeln(Dat,'Zonas de interés: ',NZonas:2);
    Writeln(Dat,'            ',Canali:6,CanalF:6);
    Writeln(Dat);
    Writeln(Dat,'Parámetros de calibración                      Código      Lambda');
    Writeln(Dat,' Ganancia:     ',P0[1]:13:11,' ± ',PE[1]:13:11,Cod[1]:6,Lam[1]:13:2);
    Writeln(Dat,'     Cero:     ',P0[2]:14:11,' ± ',PE[2]:13:11,Cod[2]:6,Lam[2]:13:2);
    Writeln(Dat);
    Writeln(Dat,'Espesores del detector                         Código      Lambda');
    Writeln(Dat,'Be:',P0[3]:27:11,' ± ',PE[3]:13:11,Cod[3]:5,Lam[3]:13:3);
    Writeln(Dat,'Ni:',P0[4]:27:13,' ± ',PE[4]:15:13,Cod[4]:3,Lam[4]:13:3);
    Writeln(Dat,'Si(Capa muerta): ',P0[5]:14:12,' ± ',PE[5]:14:12,Cod[5]:4,Lam[5]:13:3);
    Writeln(Dat);
    Writeln(Dat,'Escala de fondo                                Código      Lambda');
    Writeln(Dat,P0[6]:15:9,' ± ',PE[6]:11:9,Cod[6]:22,Lam[6]:13:2);
    Writeln(Dat,P0[7]:15:9,' ± ',PE[7]:11:9,Cod[7]:22,Lam[7]:13:2);
    Writeln(Dat,P0[8]:15:9,' ± ',PE[8]:11:9,Cod[8]:22,Lam[8]:13:2);
    Writeln(Dat);
    Writeln(Dat,'Respuesta del detector                         Código      Lambda');
    Writeln(Dat,' sigma(E)=Sqrt(Noise^2+Fano*epsilon*E)');
    Writeln(Dat,'Noise:',P0[9]:23:12,' ± ',PE[9]:14:12,Cod[9]:3,Lam[9]:13:2);
    Writeln(Dat,'Factor de Fano:',P0[10]:15:12,' ± ',PE[10]:14:12,Cod[10]:3,Lam[10]:13:2);
    Writeln(Dat);
    Writeln(Dat,'Picos detectados:       E             Fracción                   Escalon               Beta           Borde  Descripción   Códigos   Lambdas');
    l:=NPG;
    For i:=1 To NDet Do Begin
        Writeln(Dat,Sim[Zdet[i]],':',Nlin[i]:2);
        For j:=1 To Nlin[i] Do Begin
            Inc(l);
            Writeln(Dat,Perf[i,j],P0[l]:23:5,' ± ',PE[l]:7:5,P0[NL+l]:10:7,' ± ',PE[NL+l]:9:7,P0[2*NL+l]:10:6,' ± ',PE[2*NL+l]:9:7,
            P0[3*NL+l]:10:6,' ± ',PE[3*NL+l]:9:7,'   ',Borde[i,j]:2,Desc[l-NPG]:7,Cod[l]:5,Cod[l+NL]:2,Cod[l+2*NL]:2,Cod[l+3*NL]:2,Lam[l]:9:2)
        End;
    Writeln(Dat);
    End;
    Writeln(Dat,'Grupos de picos         Factores de escala             Código   Lambda');
    For i:=1 To NG Do Writeln(Dat,Grupo[i],P0[NPG+4*NL+i]:32:8,' ± ',PE[NPG+4*NL+i]:3:8,Cod[NPG+4*NL+i]:8,Lam[NPG+4*NL+i]:10:2);
    Writeln(Dat);
    Writeln(Dat,'Z     Concentraciones estimadas  Código    Lambda  (Primero los detectados)');
    For i:=NPG+NG+4*NL+1 To NPG+NG+4*NL+Ndet Do Writeln(Dat,Zdet[i-NPG-NG-4*NL]:2,P0[I]:15:8,' ± ',PE[i]:10:8,Cod[i]:6 ,Lam[i]:12:2);
    For i:=NPG+NG+4*NL+Ndet+1 To NPG+NG+4*NL+Nel Do Writeln(Dat,Znodet[i-NPG-NG-4*NL-Ndet]:2,P0[I]:15:8,' ± ',PE[i]:10:8,Cod[i]:6 ,Lam[i]:12:2);
    WriteLn(Dat);
    WriteLn(Dat,'Tolerancia para las iteraciones: ',Titer);
    Writeln(Dat);
    Writeln(Dat, 'Picos suma 	        		          Código     Lambda');
    Writeln(Dat,P0[11]:15:9,' ± ',PE[11]:11:9,Cod[11]:22,Lam[11]:13:2);
    Writeln(Dat,P0[12]:15:9,' ± ',PE[12]:11:9,Cod[12]:22,Lam[12]:13:2);
    Writeln(Dat,P0[13]:15:9,' ± ',PE[13]:11:9,Cod[13]:22,Lam[13]:13:2);
    Writeln(Dat);
    Writeln(Dat, 'Pico de Fluorescencia interna                   Código     Lambda');
    Writeln(Dat,P0[14]:15:9,' ± ',PE[14]:11:9,Cod[14]:22,Lam[14]:13:2);
    Writeln(Dat);
    WriteLn(Dat,'Espesor de la muestra (g/cm^2): ',Espesor);
    Writeln(Dat);

    {Aca empezamos a escribir la parte del film}
    If Nfilm >0 Then Begin
        Writeln(Dat,'Número total de elementos del film: ',Nfilm:2);
        Writeln(Dat,'Número de elementos detectables en el film: ',Nfdet:2);
	Write(Dat,'Números atómicos de los detectables del film: ');
	For i:=Ndet+1 To Ndet+Nfdet Do Write(Dat,Zdet[i]:3);
	Writeln(Dat);
	Write(Dat,'Números atómicos de los no etectables del film: ');
	For i:=Nnodet+1 To Nnodet+Nfnodet Do Write(Dat,Znodet[i]:3);
	Writeln(Dat);
	Writeln(Dat);
	Writeln(Dat,'Picos del film:       E             Fracción                   Escalon               Beta           Borde  Descripción   Códigos   Lambdas');
	l:=NPG;
	For i:=Ndet+1 To Ndet+Nfdet Do Begin
	    Writeln(Dat,Sim[Zdet[i]],':',Nlin[i]:2);
	    For j:=1 To Nlin[i] Do Begin
		Inc(l);
		Writeln(Dat,Perf[i,j],P0[4*NL+NG+Nel+l]:23:5,' ± ',PE[4*NL+NG+Nel+l]:7:5,P0[4*NL+NG+Nel+NLF+l]:10:7,' ± ',PE[4*NL+NG+Nel+NLF+l]:9:7,
                P0[4*NL+NG+Nel+2*NLF+l]:10:7,' ± ',PE[4*NL+NG+Nel+2*NLF+l]:9:7,P0[4*NL+NG+Nel+3*NLF+l]:10:7,' ± ',PE[4*NL+NG+Nel+3*NLF+l]:9:7,
		' ',Borde[i,j]:2,Desc[NL+l-NPG]:7,Cod[4*NL+NG+Nel+l]:5,Cod[4*NL+NG+Nel+NLF+l]:2,Cod[4*NL+NG+Nel+2*NLF+l]:2,Cod[4*NL+NG+Nel+3*NLF+l]:2,
                Lam[4*NL+NG+Nel+l]:9:2)
	    End;
	    Writeln(Dat);
	End;
	Writeln(Dat,'Z     Concentraciones estimadas  Código    Lambda  (Primero los detectados)');
	For i:=NPG+NG+4*NL+Nel+4*NLF+1 To NPG+NG+4*NL+Nel+4*NLF+Nfdet Do Writeln(Dat,Zdet[i-NPG-NG-4*NL-Nel-4*NLF+Ndet]:2,P0[i]:15:8,' ± ',PE[i]:10:8,Cod[i]:6 ,Lam[i]:12:2);
	For i:=NPG+NG+4*NL+Nel+4*NLF+Nfdet+1 To NPG+NG+4*NL+Nel+4*NLF+Nfilm Do Writeln(Dat,Znodet[i-NPG-NG-4*NL-Nel-4*NLF-Nfdet+Nnodet]:2,P0[i]:15:8,' ± ',PE[i]:10:8,Cod[i]:6 ,Lam[i]:12:2);
	WriteLn(Dat);
	WriteLn(Dat,'Espesor del film (g/cm^2): ',Espesorfilm);
    End;
    Close(Dat);
			
End;


Procedure Graficar;
Var
    file1         :text;
    i  		  :Integer;
    DiPlot        :String;
Begin
    Case SistOp of
         1: DiPlot := DiUser + 'pampa/';
         2: DiPlot := DiUser + 'pampa\';
    End;
    Assign (file1, DiPlot+'espectro.plot');
    Rewrite(file1);
    WriteLn (file1, 'set xlabel "Energy [keV]"');
    WriteLn (file1, 'set ylabel "Spectrum [number of counts]"');
    WriteLn (file1, 'set autoscale');
    WriteLn (file1, 'set log y');
    WriteLn (file1, 'set yrange [1:]');
    WriteLn (file1, 'set grid 					# para habilitar la grid (grilla)');
    WriteLn (file1, 'set xzeroaxis lt -1 linecolor rgb "black"	# para dibujar el eje x en y=0');
    WriteLn (file1, 'set key box linestyle 1 linecolor rgb "black"	# para dibujar un recuadro en la leyenda');
    WriteLn (file1, '');
    WriteLn (file1, '# para dibujar el residuo movido en un 10% del valor m\E1ximo de las cuentas del espectro:');
    WriteLn (file1, 'max(a,b)=a>b?a:b ');
    WriteLn (file1, 'a=0.');
    WriteLn (file1, '#');
    WriteLn (file1, '');
    {comienza el plot de los datos en espectro.sal}
    WriteLn (file1, '###');
    WriteLn (file1, 'plot ''' + DiSal+'.sal'' ' + 'using 1:2 title ''Spectrum'' with points pointtype 7 pointsize 0.25 linecolor rgb "black", \');
    WriteLn (file1, '     ''' + DiSal+'.sal'' ' + 'using 1:3 title ''PAMPA'' with lines linecolor rgb "black", \');
    WriteLn (file1, '     ''' + DiSal+'.sal'' ' + 'using 1:' ,4+1, ' title ''Lines'' with lines linecolor rgb "blue", \');
    For i:=2 To NL+NLF Do WriteLn (file1, '     ''' + DiSal+'.sal'' ' + 'using 1:' ,4+i, ' notitle ''Linea ',Desc[i],''' with lines linecolor rgb "blue", \'); {el 3 hace esquivar las 3 primeras columnas, que no son las líneas}
    WriteLn (file1, '     ''' + DiSal+'.sal'' ' + 'using 1:' ,4+NL+NLF+1,' title ''Asymmetries'' with lines linecolor rgb "cyan" , \');
    For i:=NL+NLF+2 To 2*NL+2*NLF+1 Do WriteLn (file1, '     ''' + DiSal+'.sal'' ' + 'using 1:' ,4+i,' with lines linecolor rgb "cyan" notitle , \');
    if (Npicosuma>0) then WriteLn (file1, '     ''' + DiSal+'.sal'' ' + 'using 1:' ,4+2*NL+NLF+2,' title ''Sum peaks'' with lines linecolor rgb "green" , \');
    For i:=2*NL+NLF+2 To 2*NL+NLF+1+Npicosuma Do WriteLn (file1, '     ''' + DiSal+'.sal'' ' + 'using 1:' ,4+i,' with lines linecolor rgb "green" notitle , \');
    WriteLn (file1, '     ''' + DiSal+'.sal'' ' + 'using 1:4 title ''Background'' with lines linecolor rgb "magenta", \');
    WriteLn (file1, '     ''' + DiSal+'.sal'' ' + 'using 1:(a=max(a,$2),0/0) notitle, '''' using 1:($3-$2-a/10) title ''Residue'' with linespoints pointsize 0.1 linecolor rgb "red" , \');
    WriteLn (file1, '     0.-a/10 with lines linetype 7 linecolor rgb "red" notitle'); {es el cero del residuo}

    {escribe las líneas necesarias para dibujar flechas sin cabeza verticales para cada línea encontrada}
    WriteLn (file1);
    WriteLn (file1, '#flechas sin cabeza donde est\E1n las lineas encontradas');

    WriteLn (file1, 'replot');
    WriteLn (file1, 'pause-1');
    close (file1);
    {grafica el archivo Ndat.sal}
    gnuplotlin;

End;


Procedure Archivos;
    Var 
    T           :Text;
    i,j,O       :Integer;
Begin
     Case SistOp of
          1: Pref := DiUser+'dat/';
          2: Pref := DiUser+'Dat\';
     End;
    Assign(T,Pref+'SIMB.ASC'); Reset(T);
    For i:=1 To 92 Do ReadLn(T,Sim[i]);
    Close(T);
    Assign(T,Pref+'OM3.ASC'); Reset(T);
    For i:=1 To 92 Do Begin
        For j:=1 To 9 Do Read(T,Omega[i,j]);
        Readln(T)
    End;
    Close(T);
    Assign(T,Pref+'BORDES.ASC');Reset(T);
    For i:=1 To 92 Do Begin
	ReadLn(T,O);
	ReadLn(T,Ek[i]);
	ReadLn(T,El1[i],El2[i],El3[i]);
	ReadLn(T,Em1[i],Em2[i],Em3[i],Em4[i],Em5[i])
    End;
    Close(T);
    Assign(T,Pref+'AES.ASC');Reset(T);  {leemos la masa atomica (g/mol)}
    For i:=1 To 92 Do  ReadLn(T,A[i]);
    Close(T);
    Assign(T,Pref+'RHO.ASC'); Reset(T); {leemos la densidad (g/cm^3)}
    For i:=1 To 92 Do ReadLn(T,Rho[i]);
    Assign(T,Pref+'CK.ASC'); Reset(T); {leemos las probabilidades de transicion de las Coster Kroning)}
    For i:=1 To 92 Do ReadLn(T,CK[i,1],CK[i,2],CK[i,3]); {f12,f13,f23}
    Close(T);
    Assign(T,Pref+'Efi.dat');Reset(T);
    For i:=1 To  26 Do Begin
         ReadLn(T, Efic[i,1],Efic[i,2]);
    End;
    Close(T);
End;


Function Func(P :DoubleArrayNP; NPR :Double) :Double;
    Var
        i,u,j,l,m               :Integer;
        En,dEsuma,Eesc,Esi      :Double;
        Sum,Fat     	        :Double;
        min,Sp,depth            :Double;
Begin
     {Creamos el arreglo X donde juntamos los parametros de P0 y los parametros refinados P}
     For i:=1 To NPG+4*NL+NG+Nel+4*NLF+Nfilm Do Begin
         If Sumcod[I]>0 Then X[i]:=(1.0-Cod[i])*P0[i]+Cod[i]*P[Sumcod[i]]
         Else X[i]:=P0[i]
     End;
     {Guardamos las nuevas concentraciones}
     For i:=1 To Nel Do C[i]:=X[NPG+4*NL+NG+i];
     l:=0;
     For i:=Nel+1 To Nel+Nfilm Do Begin
         Inc(l);
         C[i]:=X[NPG+4*NL+NG+Nel+4*NLF+l];
     End;
     {Guardamos las nuevas energias y fracciones de linea}
     l:=0;
     For i:=1 To Ndet Do For j:=1 To Nlin[i] Do Begin
         Inc(l);
         E[i,j]:=X[NPG+l];
	 Fr[i,j]:=X[NPG+NL+l];
         Tail[i,j]:=X[NPG+2*NL+l];
         Beta[i,j]:=X[NPG+3*NL+l];
         Sigma2[i,j]:=X[9]*X[9]+0.00376*X[10]*E[i,j];
     End;
     l:=0;
     For i:=Ndet+1 To Ndet+Nfdet Do For j:=1 To Nlin[i] Do Begin
         Inc(l);
         E[i,j]:=X[NPG+4*NL+NG+Nel+l];
         Fr[i,j]:=X[NPG+4*NL+NG+Nel+NLF+l];
         Tail[i,j]:=X[NPG+4*NL+NG+Nel+2*NLF+l];
         Beta[i,j]:=X[NPG+4*NL+NG+Nel+3*NLF+l];
         Sigma2[i,j]:=X[9]*X[9]+0.00376*X[10]*E[i,j];
     End;
     {Guardamos el nuevo factor de escala}
     If(ff=False) Then For i:=1 To NG Do Fesc[i]:=X[NPG+4*NL+i] {DESCOMENTAR PARA AJUSTAR FACTORES DE ESCALA DIFERENTES}
     Else If (ff=True) Then Begin
          For i:=1 To NG-1 Do Fesc[i]:=X[NPG+4*NL+1];
          Fesc[5]:=X[NPG+4*NL+5];
     End;
     {Guardamos el stopping power promedio de la muestra}
    If NFilm=0 Then Begin
         For i:=1 To Nint Do Begin
             Stopp[i]:=SPM(Eint[i])*1000;
         End;
       depth:=rhoZmax();
    End;
     {Comenzamos a evaluar la funcion intensidad}
     Sum:=0.0;
     Picosuma;
     If NFilm=0 Then Begin
        For u:=1 To Ndet Do For j:=1 To Nlin[u] Do Intens[u,j]:=Intensidad_mg(u,j,E[u,j]);
     End;
     For i:=Canali To Canalf Do Ycal[i]:=0.0;
     For i:=Canali To Canalf Do Begin
         En:=(X[1]*i+X[2]);
         If (Detector<5) Then Begin
            For u:=1 To Ndet Do For j:=1 To Nlin[u] Do Begin  {Primero sumamos la intensidad de cada linea para el canal I}
	        If NFilm=0 Then Begin  {Intensidad de las lineas cuando no hay film}
                     if Abs(E[u,j]-En)<0.3 Then Begin
                        Sp:=Perfil_Gaus(En,E[u,j],Sigma2[u,j])/(1+Tail[U,J]*Beta[U,J]);
                        Ylinea[i,u,j]:=Intens[u,j]*(1.0-X[11]*Suma[U,J])*Sp;
                     End
                     Else Begin
                          Ylinea[i,u,j]:=0.0;
                     End;
                     Sp:=Tail[U,J]*Cola(En,E[U,J],Sqrt(sigma2[u,j]),Beta[U,J])/(1+Tail[U,J]*Beta[U,J]);
                     Scola[i,U,J]:=Sp;
                     Ycola[i,u,j]:=Intens[u,j]*(1.0-X[11]*Suma[U,J])*Sp*(1.0-Resc(u,j));
		     Ycal[i]:= Ycal[i]+Ylinea[i,u,j]+Ycola[i,u,j];
	        End
		Else Begin  {Intensidad de las lineas del sustrato cuando hay film (se agrega la atenuacion de los fotones en el film)}
	             Fat:=1.0; {Factor de atenuacion}
		     For l:=Nel+1 To Nel+Nfilm Do Fat:=Fat*Exp(-1.0*C[L]*Mu(Z[L],E[U,J])/Sin((90-Toff)/180.0*Pi)*Espesorfilm);
		     Ylinea[i,u,j]:=Intensidad(u,j,E[u,j],S[u,j])*(1.0-X[11]*Suma[U,J])*
		     Perfil_Gaus(En,E[u,j],sigma2[u,j])*Fat;
		     Ycal[i]:= Ycal[i]+Ylinea[i,u,j];
		End;
            End;
            If NFilm<>0 Then Begin
               For u:=Ndet+1 To Ndet+Nfdet Do For j:=1 To Nlin[u] Do Begin {sumamos la intensidad de las lineas del film}
                     if Abs(E[u,j]-En)<0.3 Then Begin
                        Sp:=Perfil_Gaus(En,E[u,j],Sigma2[u,j])/(1+Tail[U,J]*Beta[U,J]);
                        Ylinea[i,u,j]:=Intensidad_film(u,j,E[u,j],S[u,j])*(1.0-X[11]*Suma[U,J])*Sp*(1.0-Resc(u,j));
                     End
                     Else Begin
                          Ylinea[i,u,j]:=0.0;
                     End;
                     Sp:=Tail[U,J]*Cola(En,E[U,J],Sqrt(sigma2[u,j]),Beta[U,J])/(1+Tail[U,J]*Beta[U,J]);
                     Scola[i,U,J]:=Sp;
                     Ycola[i,u,j]:=Intensidad_film(u,j,E[u,j],S[u,j])*(1.0-X[11]*Suma[U,J])*Sp*(1.0-Resc(u,j));
		     Ycal[i]:= Ycal[i]+Ylinea[i,u,j]+Ycola[i,u,j];
               End;
            End;
            If NFilm=0 then Yfondo[i]:=FondoProm(I,depth)*Efi(En)
            Else Yfondo[i]:=Fondo(En,X[6],X[7],X[8])*Fat*Efi(En);
            Ycal[i]:=Ycal[i]+Yfondo[i];  {Sumamos el fondo en el canal I}
            For m:=1 To Npicosuma Do Begin
		       dEsuma:=Esuma[m];
		       Ypicosuma[i,m]:=(ISuma[m]*X[11])*Perfil_RAS(En,dEsuma,Sigmasuma2[m]);
		       Ycal[i]:=Ycal[i]+Ypicosuma[i,m]; {Le agregamos los picos suma}
	    End;
	    If NFilm<>0 Then Begin
                For u:=Ndet+1 To Ndet+Nfdet Do For j:=1 To Nlin[u] Do Begin
		    If (E[u,j]>1.84) Then Begin
	               Eesc:=E[u,j]-1.74;  	{le agregamos los picos de escape}
		       Ycal[i]:=Ycal[i]+Intensidad_film(u,j,E[u,j],S[u,j])*Resc(u,j)*Perfil_Gaus(En,Eesc,X[9]*X[9]+0.00376*X[10]*Eesc);
		    End;
	        End;
	    End
            Else If NFilm=0 Then Begin
                 For u:=1 To Ndet Do For j:=1 To Nlin[u] Do Begin
		    If (E[u,j]>1.84) Then Begin
	               Eesc:=E[u,j]-1.74;  	{le agregamos los picos de escape}
		       Ycal[i]:=Ycal[i]+Intens[u,j]*Resc(u,j)*Perfil_Gaus(En,Eesc,X[9]*X[9]+0.00376*X[10]*Eesc);
		    End;
	        End;
            End;
            If (En>1.0) And (En<2.5) Then Begin
               Esi:=1.739;
               Ycal[i]:=Ycal[i]+X[14]*Perfil_Gaus(En,Esi,X[9]*X[9]+0.00376*X[10]*1.739);
            End;

        End;
	If Ymed[i]>0 Then Sum:=Sum+Sqr(Ymed[i]-Ycal[i])/Ymed[i]
	Else Sum:=Sum+Sqr(Ymed[i]-Ycal[i])
     End;
     {Funcion chi^2}
     Func:=Sum/(CanalF-Canali+1-NPR);
     {Funcion penalty (para que no de valores negativos)}
     For i:=8 To NPG+4*NL+NG+Nel+4*NLF+Nfilm Do Begin
	min:=0.0;
	If X[I]<0.0 Then min:=X[i];
	Func:=Func+10e10*min*min;
     End;

End;

Procedure Errores;
Type
   pNN = ^DoubleArrayNPbyNP;
    DoubleArrayNPbyNP = ARRAY [1..Np,1..np] OF Double;
Var
  Vx     :pNN;
  Guardo :Integer;
  Corr   :Text;
  i,j    :Integer;

Procedure MatDeriv;   { Calcula las dy_i/dx_k y luego la inversa de Vx }
Var
        i,u,j,l,m,k                           :Integer;
        En,dEsuma,Eesc                        :Double;
        Sum,Fat     	                      :Double;
        min,Sp2,depth                          :Double;
        Mat                                   :^DoubleArrayNp;
        XX,Yinc,Yfond,EnE                     :Double;
        EsumaE,IsumaE,S2E                     :Array[1..300] Of Double;
        SumaE                                 :Array[1..30,1..30] Of Double;
        EE,FrE,BetaE,TailE,Sigma2E            :Array[1..100,1..30] Of Double;
        CE                                    :Array[1..100] Of Double;
        FescE                                 :Array[1..10] Of Double;
Begin
     New(Mat);
     {Creamos el arreglo X donde juntamos los parametros de P0 y los parametros refinados P}
     {For i:=1 To NPG+4*NL+NG+Nel+4*NLF+Nfilm Do Begin
         If Sumcod[I]>0 Then X[i]:=(1.0-Cod[i])*P0[i]+Cod[i]*P_final[Sumcod[i]]
         Else X[i]:=P0[i]
     End;
     {ShowMessage('Al principio de func'+FloatToStr(X[13])); }
     {Guardamos las nuevas concentraciones}
     For i:=1 To Nel Do C[i]:=X[NPG+4*NL+NG+i];
     l:=0;
     For i:=Nel+1 To Nel+Nfilm Do Begin
         Inc(l);
         C[i]:=X[NPG+4*NL+NG+Nel+4*NLF+l];
     End;
     {Guardamos las nuevas energias y fracciones de linea}
     l:=0;
     For i:=1 To Ndet Do For j:=1 To Nlin[i] Do Begin
         Inc(l);
         E[i,j]:=X[NPG+l];
	 Fr[i,j]:=X[NPG+NL+l];
         Tail[i,j]:=X[NPG+2*NL+l];
         Beta[i,j]:=X[NPG+3*NL+l];
         Sigma2[i,j]:=X[9]*X[9]+0.00376*X[10]*E[i,j];
     End;
     l:=0;
     For i:=Ndet+1 To Ndet+Nfdet Do For j:=1 To Nlin[i] Do Begin
         Inc(l);
         E[i,j]:=X[NPG+4*NL+NG+Nel+l];
         Fr[i,j]:=X[NPG+4*NL+NG+Nel+NLF+l];
         Tail[i,j]:=X[NPG+4*NL+NG+Nel+2*NLF+l];
         Beta[i,j]:=X[NPG+4*NL+NG+Nel+3*NLF+l];
         Sigma2[i,j]:=X[9]*X[9]+0.00376*X[10]*E[i,j];
     End;
     {Guardamos el nuevo factor de escala}
     For i:=1 To NG Do Fesc[i]:=X[NPG+4*NL+i]; {DESCOMENTAR PARA AJUSTAR FACTORES DE ESCALA DIFERENTES}
     {For i:=1 To NG-1 Do Fesc[i]:=X[NPG+4*NL+1];
     Fesc[5]:=X[NPG+4*NL+5]; }                                    }
     {para guardar los arreglos de informacion de los picos suma}
     If P0[11]<>0.0 Then Begin
     	Picosuma;
        For M:=1 To Npicosuma Do Begin
          ESumaE[M]:=ESuma[M];
          S2E[M]:=SigmaSuma2[M];
          ISumaE[M]:=ISuma[M]
        End;
        For M:=1 To NDet Do For J:=1 To Nlin[J] Do SumaE[M,J]:=Suma[M,J]
     End
     Else For I:=1 To Ndet Do For J:=1 To Nlin[I] Do Suma[I,J]:=0.0;

     For L:=1 to NPR Do For K:=1 To NPR Do Vx^[K,L]:=0.0;

     {*****************Calculo de la intensidad de las líneas caracteristicas*******************************}
     {Guardamos el stopping power promedio de la muestra}
     If NFilm=0 Then Begin
          For i:=1 To Nint Do Begin
              Stopp[i]:=SPM(Eint[i])*1000;
          End;
        depth:=rhoZmax();
     End;
     {Comenzamos a evaluar la funcion intensidad}
     Sum:=0.0;
     For i:=Canali To Canalf Do Begin
         En:=(P0[1]*i+P0[2]);
         For K:=1 To NPG+4*NL+NG+Nel+4*NLF+Nfilm Do If Cod[K]=1 Then Begin                          { SumCod[K] SEÑALA EL PARAMETRO }
             XX:=P0[K];         { Guardo el parametro sin incrementar }
             P0[K]:=P0[K]*1.001;
             L:=0;
             If K In [1,2] Then Begin
                EnE:=En;
                En:=P0[1]*i+P0[2]
             End;
             If K In[9..13] Then Picosuma;
             If K In[NPG..NPG+2*NL]  Then For M:=1 to Ndet do For J:=1 To Nlin[M] Do Begin
                Inc(L);
                EE[m,j]:=E[m,j];
                E[m,j]:=P0[NPG+l];
                FrE[m,j]:=Fr[m,j];
                Fr[m,j]:=P0[NPG+NL+l];
                Sigma2E[M,J]:=sigma2[M,J];
                Sigma2[i,j]:=P0[9]*P0[9]+0.00376*P0[10]*E[m,j];
             End;
             If K In[NPG+2*NL+1..NPG+4*NL]  Then For M:=1 to Ndet do For J:=1 To Nlin[M] Do Begin
                Inc(L);
                TailE[m,j]:=Tail[m,j];
                Tail[m,j]:=P0[NPG+2*NL+l];
                BetaE[m,j]:=Beta[m,j];
                Beta[m,j]:=P0[NPG+3*NL+l];
             End;
             If K In[NPG+4*NL+1..NPG+4*NL+NG]  Then For m:=1 To NG Do Begin
                FescE[m]:=Fesc[m];
                Fesc[m]:=P0[NPG+4*NL+m]
             End;
             If K In[NPG+NG+4*NL+1..NPG+NG+4*NL+Nel] Then For m:=1 To Nel Do Begin
                CE[m]:=C[m];
                C[m]:=P0[NPG+4*NL+NG+m];
             End;
             If K In[NPG+4*NL+NG+Nel+1..NPG+4*NL+NG+Nel+2*NLF] Then  For m:=Ndet+1 To Ndet+Nfdet Do For j:=1 To Nlin[m] Do Begin
                Inc(L);
                EE[m,j]:=E[m,j];
                E[m,j]:=P0[NPG+4*NL+NG+Nel+l];
                FrE[m,j]:=Fr[m,j];
                Fr[m,j]:=P0[NPG+4*NL+NG+Nel+NLF+l];
                Sigma2E[M,J]:=sigma2[M,J];
                Sigma2[i,j]:=P0[9]*P0[9]+0.00376*P0[10]*E[m,j];
             End;
             If K In[NPG+4*NL+NG+Nel+2*NLF+1..NPG+4*NL+NG+Nel+4*NLF] Then For m:=Ndet+1 To Ndet+Nfdet Do For j:=1 To Nlin[m] Do Begin
                Inc(L);
                TailE[m,j]:=Tail[m,j];
                Tail[m,j]:=P0[NPG+4*NL+NG+Nel+2*NLF+l];
                BetaE[m,j]:=Beta[m,j];
                Beta[m,j]:=P0[NPG+4*NL+NG+Nel+3*NLF+l];
             End;
             If K In[NPG+4*NL+NG+Nel+4*NLF+1..NPG+4*NL+NG+Nel+4*NLF+Nfilm] Then For m:=Nel+1 To Nel+Nfilm Do Begin
                Inc(l);
                CE[m]:=C[m];
                C[m]:=P0[NPG+4*NL+NG+Nel+4*NLF+l];
             End;
             {Acá empieza la funcion func}
             Yinc:=0.0;
             If NFilm=0 Then Begin
                For u:=1 To Ndet Do For j:=1 To Nlin[u] Do Intens[u,j]:=Intensidad_mg(u,j,E[u,j]);
             End;
             If (Detector<4) Then Begin
                For u:=1 To Ndet Do For j:=1 To Nlin[u] Do Begin  {Primero sumamos la intensidad de cada linea para el canal I}
	            If NFilm=0 Then Begin  {Intensidad de las lineas cuando no hay film}
                       if Abs(E[u,j]-En)<0.3 Then Begin
                          Sp2:=Perfil_Gaus(En,E[u,j],Sigma2[u,j])/(1+Tail[U,J]*Beta[U,J]);
                          Yinc:=Yinc+Intens[u,j]*(1.0-P0[11]*Suma[U,J])*Sp2;
                       End
                       Else Begin
                          Sp2:=0.0;
                       End;
                       Sp2:=Tail[U,J]*Cola(En,E[U,J],Sqrt(sigma2[u,j]),Beta[U,J])/(1+Tail[U,J]*Beta[U,J]);
                       Yinc:=Yinc+Intens[u,j]*(1.0-P0[11]*Suma[U,J])*Sp2*(1.0-Resc(u,j));
                    End
		    Else Begin  {Intensidad de las lineas del sustrato cuando hay film (se agrega la atenuacion de los fotones en el film)}
	               Fat:=1.0; {Factor de atenuacion}
		       For l:=Nel+1 To Nel+Nfilm Do Fat:=Fat*Exp(-1.0*C[L]*Mu(Z[L],E[U,J])/Sin((90-Toff)/180.0*Pi)*Espesorfilm);
		       Yinc:=Yinc+Intensidad(u,j,E[u,j],S[u,j])*(1.0-P0[11]*Suma[U,J])*
		       Perfil_Gaus(En,E[u,j],sigma2[u,j])*Fat;
		    End;
                End;
                If NFilm<>0 Then Begin
                   For u:=Ndet+1 To Ndet+Nfdet Do For j:=1 To Nlin[u] Do Begin {sumamos la intensidad de las lineas del film}
	               Yinc:=Yinc+Intensidad_film(u,j,E[u,j],S[u,j])*(1.0-P0[11]*Suma[U,J])*Perfil_Gaus(En,E[u,j],sigma2[u,j])*
                       (1.0-Resc(u,j));
                   End;
                End;
                If NFilm=0 then Yfond:=FondoProm(I,depth)*Efi(En)
                Else Yfond:=Fondo(En,P0[6],P0[7],P0[8])*Fat*Efi(En);
                Yinc:=Yinc+Yfond;  {Sumamos el fondo en el canal I}
                For m:=1 To Npicosuma Do Begin
		       dEsuma:=Esuma[m];
		       Ypicosuma[i,m]:=(ISuma[m]*P0[11])*Perfil_RAS(En,dEsuma,Sigmasuma2[m]);
		       Yinc:=Yinc+Ypicosuma[i,m]; {Le agregamos los picos suma}
	        End;
	        If NFilm<>0 Then Begin
                   For u:=Ndet+1 To Ndet+Nfdet Do For j:=1 To Nlin[u] Do Begin
		       If (E[u,j]>1.84) Then Begin
	                  Eesc:=E[u,j]-1.74;  	{le agregamos los picos de escape}
		          Yinc:=Yinc+Intensidad_film(u,j,E[u,j],S[u,j])*Resc(u,j)*Perfil_Gaus(En,Eesc,P0[9]*P0[9]+0.00376*P0[10]*Eesc);
		       End;
	           End;
	        End
                Else If NFilm=0 Then Begin
                   For u:=1 To Ndet Do For j:=1 To Nlin[u] Do Begin
		       If (E[u,j]>1.84) Then Begin
	                  Eesc:=E[u,j]-1.74;  	{le agregamos los picos de escape}
		          Yinc:=Yinc+Intens[u,j]*Resc(u,j)*Perfil_Gaus(En,Eesc,P0[9]*P0[9]+0.00376*P0[10]*Eesc);
		       End;
	           End;
                End;
             End;
             {ShowMessage(FloatToStr(Yinc));}

             {calculo de la matriz de corelacion}

              Mat^[SumCod[K]]:=(Yinc-Ycal[i])*1000.0/XX;   {M_k=dy_i/dx_k}
              P0[K]:=XX;       {Restituyo el valor que tenia antes de derivar}
              If K In [1,2] Then En:=EnE;
              If K In[NPG..NPG+2*NL]  Then For M:=1 to Ndet do For J:=1 To Nlin[M] Do Begin
                 E[m,j]:=EE[m,j];
                 Fr[m,j]:=FrE[m,j];
                 Sigma2[m,j]:=sigma2E[m,j];
              End;
              If K In[NPG+2*NL+1..NPG+4*NL]  Then For M:=1 to Ndet do For J:=1 To Nlin[M] Do Begin
                 Tail[m,j]:=TailE[m,j];
                 Beta[m,j]:=BetaE[m,j];
              End;
              If K In[NPG+4*NL+1..NPG+4*NL+NG]  Then For m:=1 To NG Do Fesc[m]:=FescE[m];
              If K In[NPG+NG+4*NL+1..NPG+NG+4*NL+Nel] Then For m:=1 To Nel Do C[m]:=CE[m];
              If K In[NPG+4*NL+NG+Nel+1..NPG+4*NL+NG+Nel+2*NLF] Then  For m:=Ndet+1 To Ndet+Nfdet Do For j:=1 To Nlin[m] Do Begin
                 E[m,j]:=EE[m,j];
                 Fr[m,j]:=FrE[m,j];
                 Sigma2[m,j]:=sigma2E[m,j];
              End;
              If K In[NPG+4*NL+NG+Nel+2*NLF+1..NPG+4*NL+NG+Nel+4*NLF] Then For m:=Ndet+1 To Ndet+Nfdet Do For j:=1 To Nlin[m] Do Begin
                 Tail[m,j]:=TailE[m,j];
                 Beta[m,j]:=BetaE[m,j];
              End;
              If K In[NPG+4*NL+NG+Nel+4*NLF+1..NPG+4*NL+NG+Nel+4*NLF+Nfilm] Then For m:=Nel+1 To Nel+Nfilm Do C[m]:=CE[m];
              If P0[11]<>0.0 Then Begin
                 For m:=1 To Npicosuma Do Begin
                     ESuma[m]:=ESumaE[m];
                     SigmaSuma2[m]:=S2E[m];
                     ISuma[m]:=ISumaE[m]
                 End;
                 For M:=1 To NDet Do For J:=1 To Nlin[J] Do Suma[m,j]:=SumaE[m,j]
              End;
         End;
         For L:=1 to NPR Do For K:=1 To L Do   {Vx es sim\82trica}
         If Ymed[i]<>0 Then Vx^[K,L]:=Vx^[K,L]+Mat^[K]*Mat^[L]/ymed[i] {Vx[K,L]=Sum_i(dy_i/dx_k)(dy_i/dx_L)/y_i}
     End;
     For L:=1 to NPR-1 Do For K:=L+1 To NPR Do  Vx^[K,L]:=Vx^[L,K];
     Dispose(Mat)
   End;


   Procedure Gaussj(Var a: pNN;            {Matriz a invertir}
            n: integer);       {Dimension REAL de la matriz }
   Type
       iap = ^IntegerArrayNP;
       IntegerArrayNP = ARRAY [1..np] OF integer;
   Var
      big,dum,pivinv:       double;
      i,icol,irow,j,k,l,ll: integer;
      indxc,indxr,ipiv:     iap;
   Begin
        new(indxc);
        new(indxr);
        new(ipiv);
        FOR j := 1 TO n DO ipiv^[j] := 0;
        FOR i := 1 TO n DO BEGIN
            big := 0.0;
            FOR j := 1 TO n DO
            IF ipiv^[j] <> 1 THEN
            FOR k := 1 TO n DO
            IF ipiv^[k] = 0 THEN
            IF abs(a^[j,k]) >= big THEN BEGIN
                big := abs(a^[j,k]);
                irow := j;
                icol := k
            END
            ELSE IF ipiv^[k] > 1 THEN BEGIN
                writeln('pause 1 in GAUSSJ - singular matrix');
                readln
            END;
            ipiv^[icol] := ipiv^[icol]+1;
            IF irow <> icol THEN BEGIN
               FOR l := 1 TO n DO BEGIN
                   dum := a^[irow,l];
                   a^[irow,l] := a^[icol,l];
                   a^[icol,l] := dum
               END;
            END;
            indxr^[i] := irow;
            indxc^[i] := icol;
            IF a^[icol,icol] = 0.0 THEN BEGIN
               writeln('pause 2 in GAUSSJ - singular matrix');
               readln
            END;
            pivinv := 1.0/a^[icol,icol];
            a^[icol,icol] := 1.0;
            FOR l := 1 TO n DO
            a^[icol,l] := a^[icol,l] *pivinv;
            FOR ll := 1 TO n DO
            IF ll <> icol THEN BEGIN
               dum := a^[ll,icol];
               a^[ll,icol] := 0.0;
               FOR l := 1 TO n DO
                   a^[ll,l] := a^[ll,l]-a^[icol,l] *dum;
            END
        END;
        FOR l := n DOWNTO 1 DO
        IF indxr^[l] <> indxc^[l] THEN
        FOR k := 1 TO n DO BEGIN
            dum := a^[k,indxr^[l]];
            a^[k,indxr^[l]] := a^[k,indxc^[l]];
            a^[k,indxc^[l]] := dum
        END;
        dispose(ipiv);
        dispose(indxr);
        dispose(indxc)
   End;

Begin  {Errores}
   New(Vx);
   Guardo:=Abaco;
   Abaco:=1;
   MatDeriv;         { --> Devuelve Vx^(-1) }
   GaussJ(Vx,NPR);   { Invierte Vx^(-1) }
   Abaco:=Guardo;
   Assign(Corr,DiUser+'/pampa/corr.dat');
   Rewrite(Corr);
   For i:=1 to NPR Do Begin
       For j:=1 to NPR Do
           if Vx^[i,i]*Vx^[j,j]>0 then Write(Corr,Vx^[i,j]/Sqrt(Vx^[i,i]*Vx^[j,j]):5:2,' ')
           else write('   error');
           Writeln(Corr)
   End;   {Ahora S\D6 escribe la matriz de correlaci\A2n}
   For i:=1 To NPR Do PE_ref[i]:=Sqrt(ABS(Vx^[i,i])*chi2); {El ABS deber\A1a sobrar}
   Dispose(Vx);
   Close(Corr)
End;

Function Perfil_Gaus(E: Double;Var En: Double;sigma2: Double) :Double;
Begin

    Perfil_Gaus:=(1.0/(sqrt(2*PI*sigma2)))*exp((-1.0*(E-En)*(E-En))/(2.0*sigma2));

End;


Function Perfil_RAS(E: Double;Var En: Double;sigma2: Double) :Double;
	Var
	w,d				:Double;														
Begin

	d:=X[12];
	w:=X[13];
	Perfil_RAS:=exp(d*(E-En))*(1.0/(exp((E-En)/w)+1.0));
End;


Function Intensidad(U: Integer;J: Integer; E: Double; S: Double) :Double; {(elemento,linea,energia,seccion eficaz)}
    Var    
    b          :Integer;
    alpha      :Double;
    ch1,ch2    :Char;
Begin
 
    ch1:=Borde[U,J][1];
    ch2:=Borde[U,J][2];
    {Vemos el borde de la linea para saber donde evaluar el Omega}
    If ch1='K' Then b:=1
    Else If ch1='L' Then  Begin
   		If ch2='3' Then b:=4
    	        Else If ch2='2' Then b:=3
    	        Else If ch2='1' Then b:=2;
    End
    Else If ch1='M' Then  b:=5;
    {Vemos el borde de la linea para el factor de escala}
    alpha:=0.0;
    If ch1='K' Then alpha:=Fesc[1]
    Else If ch1='L' Then Begin
		If ch2='3' Then alpha:=Fesc[2]
		Else If ch2='2' Then alpha:=Fesc[3]
		Else If ch2='1' Then alpha:=Fesc[4];
	End
    Else If ch1='M' Then  Begin
 		If ch2='5' Then alpha:=Fesc[5]
		Else If ch2='4' Then alpha:=Fesc[6]
		Else If ch2='3' Then alpha:=Fesc[7]
		Else If ch2='2' Then alpha:=Fesc[8]
		Else If ch2='1' Then alpha:=Fesc[9]
                Else If ch2='M' Then alpha:=Fesc[5];
    End;
    Intensidad:=0.602214e23*Fr[U,J]*Omega[Zdet[U],b]*alpha*C[U]*Efi(E)*S*Espesor/A[Zdet[U]]	;
   
End;

Function Intensidad_mg(U: Integer;J: Integer; E: Double) :Double; {(elemento,linea,energia)}
    Var
    b,i          :Integer;
    alpha,h,integral_trap  :Double;
    ch1,ch2    :Char;
Begin

    ch1:=Borde[U,J][1];
    ch2:=Borde[U,J][2];
    {Vemos el borde de la linea para saber donde evaluar el Omega}
    If ch1='K' Then b:=1
    Else If ch1='L' Then  Begin
   		If ch2='3' Then b:=4
    	        Else If ch2='2' Then b:=3
    	        Else If ch2='1' Then b:=2;
    End
    Else If ch1='M' Then  b:=5;
    {Vemos el borde de la linea para el factor de escala}
     alpha:=0.0;
     If ch1='K' Then alpha:=Fesc[1]
     Else If ch1='L' Then Begin
 		If ch2='3' Then alpha:=Fesc[2]
 		Else If ch2='2' Then alpha:=Fesc[3]
 		Else If ch2='1' Then alpha:=Fesc[4];
 	End
     Else If ch1='M' Then  Begin
  		If ch2='5' Then alpha:=Fesc[5]
 		Else If ch2='4' Then alpha:=Fesc[6]
 		Else If ch2='3' Then alpha:=Fesc[7]
 		Else If ch2='2' Then alpha:=Fesc[8]
 		Else If ch2='1' Then alpha:=Fesc[9]
                Else If ch2='M' Then alpha:=Fesc[5];
     End;
    Intensidad_mg:=0.602214e23*Fr[U,J]*Omega[Zdet[U],b]*alpha*C[U]*Efi(E)/A[Zdet[U]];
    {Ahora hay que agregar la integral a intensidad_mg}
    h:=Einc/Nint;
    integral_trap:=0.5*h*(CS[Nint,u,j]*Trans(Nint,E,u,j)/Stopp[Nint]+CS[1,u,j]*Trans(1,E,u,j)/Stopp[1]);
    For i:=2 To Nint-1 Do Begin
        integral_trap:=integral_trap+h*CS[i,U,J]*Trans(i,E,u,j)/Stopp[i];
    End;
    Intensidad_mg:=Intensidad_mg*integral_trap;
End;

Function Intensidad_mg2(U: Integer;J: Integer; E: Double) :Double; {(elemento,linea,energia,seccion eficaz)}
    Var
    b,i          :Integer;
    alpha,h,integral_trap  :Double;
    ch1,ch2    :Char;
Begin

    ch1:=Borde[U,J][1];
    ch2:=Borde[U,J][2];
    {Vemos el borde de la linea para saber donde evaluar el Omega}
    If ch1='K' Then b:=1
    Else If ch1='L' Then  Begin
   		If ch2='3' Then b:=4
    	        Else If ch2='2' Then b:=3
    	        Else If ch2='1' Then b:=2;
    End
    Else If ch1='M' Then  b:=5;
    {Vemos el borde de la linea para el factor de escala}
     alpha:=0.0;
     If ch1='K' Then alpha:=Fesc[1]
     Else If ch1='L' Then Begin
 		If ch2='3' Then alpha:=Fesc[2]
 		Else If ch2='2' Then alpha:=Fesc[3]
 		Else If ch2='1' Then alpha:=Fesc[4];
 	End
     Else If ch1='M' Then  Begin
  		If ch2='5' Then alpha:=Fesc[5]
 		Else If ch2='4' Then alpha:=Fesc[6]
 		Else If ch2='3' Then alpha:=Fesc[7]
 		Else If ch2='2' Then alpha:=Fesc[8]
 		Else If ch2='1' Then alpha:=Fesc[9]
                Else If ch2='M' Then alpha:=Fesc[5];
     End;
    Intensidad_mg2:=0.602214e23*Fr[U,J]*Omega[Zdet[U],b]*alpha*C[U]*Efi(E)/A[Zdet[U]];
    h:=2*Einc/Nint;
    integral_trap:=0.5*h*(CS[Nint,u,j]*Trans2(Nint,E,u,j)/Stopp[Nint]+CS[1,u,j]*Trans2(1,E,u,j)/Stopp[1]);
    For i:=2 To Round(Nint/2)-1 Do Begin
        integral_trap:=integral_trap+h*CS[2*i,U,J]*Trans2(i,E,u,j)/Stopp[2*i];
    End;
    Intensidad_mg2:=Intensidad_mg2*integral_trap;
End;


Function Intensidad_film(U: Integer;J: Integer; E: Double; S: Double ) :Double; {(elemento,linea,energia)}
    Var    
    b          :Integer;
    alpha      :Double;
    ch1,ch2    :Char;
Begin
 
    ch1:=Borde[U,J][1];
    ch2:=Borde[U,J][2];
    {Vemos el borde de la linea para saber donde evaluar el Omega}
    If ch1='K' Then b:=1
    Else If ch1='L' Then  Begin
       		If ch2='3' Then b:=4
		Else If ch2='2' Then b:=3
		Else If ch2='1' Then b:=2;
    End
    Else If ch1='M' Then  b:=5;
       {Vemos el borde de la linea para el factor de escala}
    alpha:=0.0;
    If ch1='K' Then alpha:=Fesc[1]
    Else If ch1='L' Then Begin
		If ch2='3' Then alpha:=Fesc[2]
		Else If ch2='2' Then alpha:=Fesc[3]
		Else If ch2='1' Then alpha:=Fesc[4];
	End
    Else If ch1='M' Then  Begin
 		If ch2='5' Then alpha:=Fesc[5]
		Else If ch2='4' Then alpha:=Fesc[6]
		Else If ch2='3' Then alpha:=Fesc[7]
		Else If ch2='2' Then alpha:=Fesc[8]
		Else If ch2='1' Then alpha:=Fesc[9]
                Else If ch2='M' Then alpha:=Fesc[5];
    End;
    Intensidad_film:=0.602214e23*Fr[U,J]*Omega[Zdet[U],b]*alpha*C[U-Ndet+Nel]*Efi(E)*S*Espesorfilm/A[Zdet[U]];

End;

Function RhoZmax() :Double;
Var
   i:Integer;
   h,z:Double;
Begin
     h:=Einc/Nint;
     z:=0.5*h*(1.0/Stopp[1]+1.0/Stopp[Nint]);   {valores de los extremos de la integral}
     For i:=2 To Nint-1 Do Begin               {Calculo de la integral por método del trapezoide}
         z:=z+h*(1.0/Stopp[i]);
     End;
    RhoZmax:=z;
End;

Function Trans(k:Integer;Ener:double;u,j:Integer) :Double; {Factor de transmision para muestras gruesas}
Var
   i:Integer;
   h,Tr:Double;
Begin
     h:=Einc/Nint;
    if(k=Nint) Then Tr:=0
    else Begin
         Tr:=0.5*h*(1.0/Stopp[k]+1.0/Stopp[Nint]);   {valores de los extremos de la integral}
         For i:=k+1 To Nint-1 Do Begin               {Calculo de la integral por método del trapezoide}
 	     Tr:=Tr+h*(1.0/Stopp[i]);
         End;
    end;
    Trans:=Exp(-MuM(Ener)*(1.0/sin((90-Toff)*Pi/180.0))*Tr);
End;

Function Trans2(k:Integer;Ener:double;u,j:Integer) :Double; {Factor de transmision para muestras gruesas}
Var
   i:Integer;
   h,Tr:Double;
Begin
     h:=2*Einc/Nint;
    if(k=Nint) Then Tr:=0
    else Begin
         Tr:=0.5*h*(1.0/Stopp[k]+1.0/Stopp[Nint]);   {valores de los extremos de la integral}
         For i:=k+1 To Round(Nint/2)-1 Do Begin               {Calculo de la integral por método del trapezoide}
 	     Tr:=Tr+h*(1.0/Stopp[2*i]);
         End;
    end;
    Trans2:=Exp(-MuM(Ener)*(1.0/sin((90-Toff)*Pi/180.0))*Tr);
End;


Function Fondo(E:Double;A:Double;B:Double;C:Double) :Double; {Fondo bremsstrahlung (Energia, escala de fondo)}
Begin
	
	Fondo:=C*Exp((A-B*E)*Ln(10));

End;

Function FondoProm(I:Integer;d:Double):Double;
Var
  Epe,GGaus2,GGaus,fond,fondi,sigmagaus,En,Enm1,DE,Naoux,Sumagausiana:Double;
  k,N:Integer;

Begin
 En:=X[1]*I+X[2];
 Enm1:=X[1]*(I-1)+X[2];
 GGaus:=X[9]+X[10]*0.00376*En;
 sigmagaus:=GGaus*3/sqrt(LN(2));
 DE:=En-Enm1;
 Naoux:=Abs(sigmagaus/DE);
 N:=Round(Naoux);
 fond:=0.0;
 Sumagausiana:=0.0;
 For K:=1 To 2*N+1 Do
  Begin
    Epe:=X[1]*(I-K+N+1)+X[2];{energia donde evaluo la absocion y la gausiana};
    GGaus2:=X[9]+X[10]*0.00376*Epe;{semi ancho energetico por dispersion angular}
    if Epe<=0.1 Then fondi:=1.0
    else fondi:=FondAt(Epe,d)*Exp(-LN(2)*Sqr(En-Epe)/sqr(GGaus2))/GGaus2;
    fond:=fond+Fondo(Epe,X[6],X[7],X[8])*fondi;
    Sumagausiana:=Sumagausiana+fondi
  End;
 Fond:=Fond/Sumagausiana;
 FondoProm:=Fond*sqrt(Pi)/2.0;
 End;


Function FondAt(En:Double;depth:Double) :Double;    {atenuacion del fondo para muestras gruesas}
Begin
    FondAt:=(1.0-exp(-MuM(En)*depth/sin((90-toff)/180.0*Pi)))/(MuM2(En)/sin((90-toff)/180.0*Pi))
end;


Procedure TablaEfi();
 Var
    file1  :Text;
    I      :Integer;
    Ei     :Double;
    DiEfi  :String;
Begin
    Case SistOp of
        1: DiEfi := DiUser+ '/pampa/';
        2: DiEfi := DiUser+ '\pampa\';
    End;
    Assign (file1, DiUser+'eficiencia.dat');
    Rewrite(file1);
    Ei:=0.1;
    For I:=1 To 3000 Do Begin
        WriteLn (file1, Ei,'   ', Efi(Ei));
        Ei:=Ei+0.01;
    End;
    Close(file1);
end;



Function Efi(Ei:Double) :Double;
	Var
 	Eficien,Muwin,MuHB,E,E1,E2,ef,ef1,ef2	:Double;
         Dat             :Text;
         nl,nm,i           :Integer;
Begin

	If Detector=1 Then Begin   {Detector de Si(Li) de Porto Alegre}
            Eficien:=exp(-(Mu(4,Ei)*RBe*P0[3]+Mu(28,Ei)*RNi*P0[4]+Mu(14,Ei)*RSi*P0[5])); {Debo multiplicar el mu por la densidad}
	End;				 {mu (cm^2/g) y densidad (g/cm^3)}
	If Detector=2 Then Begin  {Detector Si(Li) de Sevilla}
            Eficien:=exp(-(Mu(4,Ei)*RBe*P0[3]+Mu(13,Ei)*RAl*P0[4]+Mu(14,Ei)*RSi*P0[5]))*(1-exp(-(Mu(14,Ei)*RSi*0.3))); {Debo multiplicar el mu por la densidad}
        End;
	If Detector=3 Then Begin {Bariloche}
	    {Muwin:=0.69*Mu(6,Ei)+0.21*Mu(8,Ei)+0.07*Mu(7,Ei)+0.03*Mu(1,Ei);
	    {MuHB:=0.085*Mu(1,Ei)+0.915*Mu(5,Ei);       }

	    {Efic:=exp(-Muwin*RMoxtek*P0[3]-MuHB*2.0*0.000002)*(0.77+0.23*exp(-Mu(14,Ei)*RSi*0.038)); }
            Efic:=exp(-Muwin*RMoxtek*P0[3])*(0.77+0.23*exp(-Mu(14,Ei)*RSi*0.038));
	    Efic:=Efic*exp(-(Mu(13,Ei)*RAl*P0[4]+Mu(14,Ei)*RSi*P0[5]))*(1-exp(-Mu(14,Ei)*RSi*0.045));   }

            {ShowMessage(FloatToStr(Ei)+'  '+FloatToStr(Efic));
     	    Assign(Dat,DiUser+'/Doctorado/programa_PAMPA/PAMPA_prueba2/dat/EficienciaBariloche.eff');
            Reset(Dat);
    	    nl:=0; {numero de lineas del archivo}
    	    nm:=0;
    	    Repeat		{vemos cuantas lineas tiene el archivo para ese elemento}
    		Readln(Dat);
    		Inc(nl);
    	    Until(SeekEOf(Dat));
    	    Reset(Dat);
    	    For I:=1 To nl Do Begin
            Readln(Dat, E,Ef);
    		If Ei>E Then Begin {aca se ve cual es la energia anterior a la einc }
    			E1:=E;
    			Ef1:=Ef;
    			Inc(nm);
    		End;
            End;
            Reset(Dat);
    	    If nm=0 Then begin  {esto es cuando la Es es menor que el primer valor de E}
    		Readln(Dat,E1,ef);
    		ef1:=ef;
    		Readln(Dat,E2,ef);
    		ef2:=ef;
            End
    	    Else If nm=nl Then begin  {esto es cuando la einc es mayor que el ultimo valor de E}
    		For I:=1 To nl-2 Do Readln(Dat);
    		Readln(Dat,E1,ef);
    		ef1:=ef;
    		Readln(Dat,E2,ef);
    		ef2:=ef;
    	    End
    	    Else Begin
    		For I:=1 To nm Do Readln(Dat);
    		Readln(Dat,E2,ef);
    		ef2:=ef;				    {se lee la energia posterior a la einc }
    	    End;
    	    Close(Dat);
    	    Efic:=((Ei-E1)*(ef2-ef1))/(E2-E1)+ef1;  }
            nm:=0;
    	    For I:=1 To 26 Do Begin
    		If Ei>Efic[i,1] Then Begin {aca se ve cual es la energia anterior a la einc }
    			E1:=Efic[i,1];
    			Ef1:=Efic[i,2];
    			Inc(nm);
    		End;
            End;
    	    If nm=0 Then begin  {esto es cuando la Es es menor que el primer valor de E}
                E1:=Efic[1,1];
    		ef1:=Efic[1,2];
                E2:=Efic[2,1];
    		ef2:=Efic[2,2];
            End
    	    Else If nm=26 Then begin  {esto es cuando la einc es mayor que el ultimo valor de E}
                E1:=Efic[25,1];
    		ef1:=Efic[25,2];
                E2:=Efic[26,1];
    		ef2:=Efic[26,2];
    	    End
            Else Begin
                E2:=Efic[nm+1,1];
    		ef2:=Efic[nm+1,2];				    {se lee la energia posterior a la einc }
    	    End;
    	    Eficien:=((Ei-E1)*(ef2-ef1))/(E2-E1)+ef1;
        End;
	If Detector=4 Then Begin {Desde archivo}
           Assign(Dat,Dieff);
           Reset(Dat);
    	   nl:=0; {numero de lineas del archivo}
    	   nm:=0;
    	   Repeat		{vemos cuantas lineas tiene el archivo}
    		Readln(Dat);
    		Inc(nl);
    	   Until(SeekEOf(Dat));
    	   Reset(Dat);
    	   For I:=1 To nl Do Begin
               Readln(Dat, E,Ef);
    	       If Ei>E Then Begin {aca se ve cual es la energia anterior a la einc }
    	          E1:=E;
    		  Ef1:=Ef;
    		  Inc(nm);
    	       End;
           End;
           Reset(Dat);
    	   If nm=0 Then begin  {esto es cuando la Es es menor que el primer valor de E}
    		Readln(Dat,E1,ef);
    		ef1:=ef;
    		Readln(Dat,E2,ef);
    		ef2:=ef;
           End
    	   Else If nm=nl Then begin  {esto es cuando la einc es mayor que el ultimo valor de E}
    		For I:=1 To nl-2 Do Readln(Dat);
    		Readln(Dat,E1,ef);
    		ef1:=ef;
    		Readln(Dat,E2,ef);
    		ef2:=ef;
    	   End
    	   Else Begin
    		For I:=1 To nm Do Readln(Dat);
    		Readln(Dat,E2,ef);
    		ef2:=ef;				    {se lee la energia posterior a la einc }
    	   End;
    	   Close(Dat);
    	   Eficien:=((Ei-E1)*(ef2-ef1))/(E2-E1)+ef1;

        end;
     Efi:=Eficien;
	
End;

Function MuM(E:Double):Double;     {Mu masico medio}
Var
     i:Integer;
Begin
     MuM:=0.0;
     For i:=1 To Nel Do Begin
         MuM:=MuM+Mu(Z[i],E)*C[i];
     End;
End;

Function MuM2(E:Double):Double;
Var
     i:Integer;
Begin
     MuM2:=0.0;
     For i:=1 To Nel Do Begin
         MuM2:=MuM2+Mu(Z[i],E)*C[i]*Rho[Z[i]];
     End;
End;


Function Mu(Z:Integer;E:Double):Double; {Heinrich} {Es el mu masico (cm^2/g)}
	VAR       
	Q,C1,C2,C3,C,N,Bias,B	:Double;
Begin
	If E>Ek[Z] Then
		If Z<6 Then Begin
			C:=-2.87536E-4 +1.808599E-3*Z;
			N:=3.34745+(0.02652873-0.01273815*Z)*Z;
			Bias:=-103.0+18.2*Z;
			B:=24.4545 +(155.6055 -14.15422*Z)*Z
		End
		Else Begin
			C:=5.253E-3+(1.33257E-3+(-7.5937E-5+(1.69357E-6-1.3975E-8*Z)*Z)*Z)*Z;
			N:=3.112-0.0121*Z;
			Bias:=0.0;
			B:=(47.0+(6.52-0.152624*Z)*Z)*Z
		End
	Else If E>El3[Z] Then Begin
		C:=-9.24E-5+(1.41478E-4+(-5.24999E-6+(9.85296E-8+(-9.07306E-10+3.19245E-12*Z)*Z)*Z)*Z)*Z;
		N:=2.7575+(1.889E-3-4.982E-5*Z)*Z;
		Bias:=0.0;
		B:=(17.8096+(6.7429E-2+(1.253775E-2-1.16286E-4*Z)*Z)*Z)*Z;
		If E<=El1[Z] Then Begin
			C:=C*0.858;
			If E<=El2[Z] Then C:=C*(0.8933+(-8.29E-3+6.38E-5*Z)*Z)
		End
	End
	Else If E>Em1[Z] Then Begin
		C:=1.889757E-2+(-1.8517159E-3+(6.9602789E-5+(-1.1641145E-6+7.2773258E-9*Z)*Z)*Z)*Z;
		N:=0.5385+(8.4597E-2+(-1.08246E-3 +4.4509E-6*Z)*Z)*Z;
		Bias:=(5.654+(-0.536839169+(1.8972278E-2-1.683474E-4*Z)*Z)*Z)*Z;
		B:=(10.2575657+(-0.822863477+(2.63199611E-2-1.8641019E-4*Z)*Z)*Z)*Z;
		If Z>29 Then Begin
			C:=3.0039E-3+(-1.73663566E-4+(4.0424792E-6+(-4.0585911E-8+1.497763E-10*Z)*Z)*Z)*Z;
			If Z>60 Then Bias:=(-1232.4022+(51.114164+(-0.699473097+3.1779619E-3*Z)*Z)*Z)*Z
		End
	End
	Else If E>Em5[Z] Then Begin
		C1:= 7.7708E-5+(-7.83544E-6+(2.209365E-7-1.29086E-9*Z)*Z)*Z;
		C2:= 1.406+(1.62E-2    +(-6.561E-4  +4.865E-6*Z)*Z)*Z;
		C3:= 0.584+(1.955E-2   -1.285E-4*Z)*Z;
		N:= 3.0-4E-3*Z;
		Bias:=(2.51+(-5.2E-2    +3.78E-4*Z)*Z)*Em4[Z]*1000;
		B:=(4.62-0.04*Z)*Z;
		If E>Em2[Z] Then C:=C1*C2*C3
    	Else If E>Em3[Z] Then C:=(1.082+1.366E-3*Z)*C1*C2
        Else If E>Em4[Z] Then C:=0.95*C1*C2
		Else C:=(1.6442+(-0.048+4.0664E-4*Z)*Z)*C1*C2
    End
    Else Begin
		C:=( 4.3156E-3 +(-1.4653E-4 +(1.707073E-6 - 6.69827E-9*Z)*Z)*Z)*1.08;
		N:=0.3736+2.401E-2*Z;
		Bias:=-113.0+4.5*Z;
		B:=(19.64+(-0.61239+5.39309E-3*Z)*Z)*Z
	End;
	Q:=SQR(1.0*Z);	
	Q:=C*EXP(N*LN(12.397/E))*SQR(Q);
	Q:=Q/A[Z];
	C:=(Bias-E*1000.0)/B;
	Mu:=-Q*(EXP(C)-1.0)
End;

Procedure TablaSP();
 Var
    T  :Text;
    I,j    :Integer;
Begin
     For i:=1 To Nel Do Begin
         Assign(T,DiUser+'/dat/SP-'+FloatToStr(Z[i])+'.dat');
         Reset(T);
         Readln(T);Readln(T);Readln(T);
         For j:=1 To  79 Do Begin
             ReadLn(T,SpAr[i,j,1],SpAr[i,j,2]);
         End;
         Close(T);
     end;
end;

Function SPM(Ei:Double)  :Double; {Ei en keV}
	Var
        Dat                       :Text;
	i,j,nm,nl	          :Integer;
	E1,E2,S1,S2,SP,E,S        :Double;
Begin
     SP:=0.0;
     For j:=1 To Nel Do Begin
	     nl:=79; {numero de lineas del archivo}
	     nm:=0;
	     For I:=1 To nl Do Begin
		      If Ei>SpAr[j,i,1] Then Begin {aca se ve cual es la energia anterior a la einc }
			     E1:=SpAr[j,i,1];
			     S1:=SpAr[j,i,2];
			     Inc(nm);
		      End;
	     End;
	     If nm=0 Then begin  {esto es cuando la Es es menor que el primer valor de E}
    	              E1:=SpAr[j,1,1];
		      S1:=SpAr[j,1,2];
    	              E2:=SpAr[j,2,1];
		      S2:=SpAr[j,2,2];
	     End
	     Else If nm=nl Then begin  {esto es cuando la einc es mayor que el ultimo valor de E}
	              E1:=SpAr[j,nl-1,1];
    	              S1:=SpAr[j,nl-1,2];
	              E2:=SpAr[j,nl,1];
    	              S2:=SpAr[j,nl,2];
	     End
	     Else Begin
     	              E2:=SpAr[j,nm+1,1];
		      S2:=SpAr[j,nm+1,2];				    {se lee la energia posterior a la einc }
	     End;
	     SP:=SP+(((Ei-E1)*(S2-S1))/(E2-E1)+S1)*C[j];

     end;
    { SP:=0.0;
     nl:=79;
     For j:=1 To Nel Do Begin
         nm:=0;
	 For I:=1 To nl Do Begin
	      {Readln(Dat, E,S); }
	      If Ei>SpAr[j,i,1] Then Begin {aca se ve cual es la energia anterior a la einc }
		     E1:=SpAr[j,i,1];
		     S1:=SpAr[j,i,2];
		     Inc(nm);
	      End;
	 End;
	 If nm=0 Then begin  {esto es cuando la Es es menor que el primer valor de E}
	    S1:=SpAr[j,1,2];
	    S2:=SpAr[j,2,2];
	End
	Else If nm=nl Then begin  {esto es cuando la einc es mayor que el ultimo valor de E}
	    S1:=SpAr[j,nl-1,2];
	    S2:=SpAr[j,nl,2];
	End
	Else Begin
             E2:=SpAr[j,nm+1,1];
	     S2:=SpAr[j,nm+1,2];				    {se lee la energia posterior a la einc }
	End;
	SP:=SP+(((Ei-E1)*(S2-S1))/(E2-E1)+S1)*C[j];
     End;           }
     SPM:=SP;
End;

Function SPM2(Ei:Double)  :Double; {Ei en keV}
	Var
	i,j,nm,nl			:Integer;
	E1,E2,S1,S2,S,E,SP      :Double;
	Dat			:Text;
Begin
     SP:=0.0;
     For j:=1 To Nel Do Begin
	     Assign(Dat,DiUser+'/dat/SP-'+FloatToStr(Z[j])+'.dat');
	     Reset(Dat);
	     nl:=0; {numero de lineas del archivo}
	     nm:=0;
	     Repeat		{vemos cuantas lineas tiene el archivo para ese elemento}
		     Readln(Dat);
		     Inc(nl);
	     Until(SeekEOf(Dat));
	     nl:=nl-3;
	     Reset(Dat);
	     Readln(Dat);Readln(Dat);Readln(Dat);
	     For I:=1 To nl Do Begin
		      Readln(Dat, E,S);
		      If Ei>E Then Begin {aca se ve cual es la energia anterior a la einc }
			     E1:=E;
			     S1:=S;
			     Inc(nm);
		      End;
	     End;
	     Reset(Dat);Readln(Dat);Readln(Dat);Readln(Dat);
	     If nm=0 Then begin  {esto es cuando la Es es menor que el primer valor de E}
		      Readln(Dat,E1,S);
		      S1:=S;
		      Readln(Dat,E2,S);
		      S2:=S;
	     End
	     Else If nm=nl Then begin  {esto es cuando la einc es mayor que el ultimo valor de E}
		      For I:=1 To nl-2 Do Readln(Dat);
		      Readln(Dat,E1,S);
		      S1:=S;
		      Readln(Dat,E2,S);
		      S2:=S;
	     End
	     Else Begin
	          For I:=1 To nm Do Readln(Dat);
		      Readln(Dat,E2,S);
		      S2:=S;				    {se lee la energia posterior a la einc }
	     End;
	     Close(Dat);
	     SP:=SP+(((Ei-E1)*(S2-S1))/(E2-E1)+S1)*C[j]*Rho[Z[j]];
     End;
     SPM2:=SP;
End;

Function SigmaK(Z:Integer;Ei:Double)  :Double; {datos de PIXE2010: Proton/alpha Ionization (k, l, m shell), Tabulated Cross Section Library.}
	Var                          {E en MeV}
	i,nm,nl			:Integer;
	E1,E2,S1,S2,E,sigma	:Double;
	Dat			:Text;
        DiSigma                 :String;
Begin
        Case SistOp of
             1: DiSigma := DiUser+'/dat/proton/k/ecpssr_he_hs_ua/k-';
             2: DiSigma := DiUser+'\dat\proton\k\ecpssr_he_hs_ua\k-';
        End;
	Assign(Dat,DiSigma+FloatToStr(Z)+'.dat');
	Reset(Dat);
	nl:=0; {numero de lineas del archivo}
	nm:=0;
	Repeat		{vemos cuantas lineas tiene el archivo para ese elemento}
		Readln(Dat);
		Inc(nl);
	Until(SeekEOf(Dat));
	nl:=nl-2; {las dos ultimas lineas no cuentan}
	Reset(Dat);
	For I:=1 To nl Do Begin
		Readln(Dat, E, sigma);
		If (Ei/1000)>E Then Begin {aca se ve cual es la energia anterior a la einc }
			E1:=E;
			s1:=sigma;
			Inc(nm);
		End;
	End;
	Reset(Dat);
	If nm=0 Then begin  {esto es cuando la einc es menor que el primer valor de E}
		Readln(Dat, E1, s1);
		Readln(Dat, E2, s2);
	End
	Else If nm=nl Then begin  {esto es cuando la einc es mayor que el ultimo valor de E}
		For I:=1 To nl-2 Do Readln(Dat);
		Readln(Dat, E1, s1);
		Readln(Dat, E2, s2);
	End
	Else Begin
		For I:=1 To nm Do Readln(Dat);
		Readln(Dat, E2, s2);    {se lee la energia posterior a la einc }
	End;
	E1:=E1*1000;
	E2:=E2*1000;

	Close(Dat);
	SigmaK:=(((Ei-E1)*(S2-S1))/(E2-E1)+S1)*1e-24;	{de mb a cm^2}

End;


Function SigmaL1(Z:Integer;Ei:Double)  :Double; {datos de PIXE2010: Proton/alpha Ionization (k, l, m shell), Tabulated Cross Section Library.}
	Var
	i,nm,nl				:Integer;
	E1,E2,S1,S2,E,sigma	        :Double;
	Dat				:Text;
        DiSigma                         :String;
Begin
        Case SistOp of
             1: DiSigma := DiUser+'/dat/proton/l/ecpssr/l1-';
             2: DiSigma := DiUser+'\dat\proton\l\ecpssr\l1-';
        End;
	Assign(Dat,DiSigma+FloatToStr(Z)+'.dat');
	Reset(Dat);
	nl:=0; {numero de lineas del archivo}
	nm:=0;
	Repeat		{vemos cuantas lineas tiene el archivo para ese elemento}
		Readln(Dat);
		Inc(nl);
	Until(SeekEOf(Dat));
	nl:=nl-2;
	Reset(Dat);
	For I:=1 To nl Do Begin
		Readln(Dat, E, sigma);
		If (Ei/1000)>E Then Begin {aca se ve cual es la energia anterior a la einc }
			E1:=E;
			s1:=sigma;
			Inc(nm);
		End;
	End;
	Reset(Dat);
	If nm=0 Then begin  {esto es cuando la einc es menor que el primer valor de E}
		Readln(Dat, E1, s1);
		Readln(Dat, E2, s2);
	End
	Else If nm=nl Then begin  {esto es cuando la einc es mayor que el ultimo valor de E}
		For I:=1 To nl-2 Do Readln(Dat);
		Readln(Dat, E1, s1);
		Readln(Dat, E2, s2);
	End
	Else Begin
		For I:=1 To nm Do Readln(Dat);
		Readln(Dat, E2, s2);    {se lee la energia posterior a la einc }
	End;
	E1:=E1*1000;
	E2:=E2*1000;

	Close(Dat);

	SigmaL1:=(((Ei-E1)*(S2-S1))/(E2-E1)+S1)*1e-24;	{de mb a cm^2}

End;

Function SigmaL2(Z:Integer;Ei:Double)  :Double; {datos de PIXE2010: Proton/alpha Ionization (k, l, m shell), Tabulated Cross Section Library.}
	Var
	i,nm,nl				:Integer;
	E1,E2,S1,S2,E,sigma	        :Double;
	Dat				:Text;
        DiSigma                         :String;
Begin
        Case SistOp of
             1: DiSigma := DiUser+'/dat/proton/l/ecpssr/l2-';
             2: DiSigma := DiUser+'\dat\proton\l\ecpssr\l2-';
        End;
	Assign(Dat,DiSigma+FloatToStr(Z)+'.dat');
	Reset(Dat);
	nl:=0; {numero de lineas del archivo}
	nm:=0;
	Repeat		{vemos cuantas lineas tiene el archivo para ese elemento}
		Readln(Dat);
		Inc(nl);
	Until(SeekEOf(Dat));
	nl:=nl-2;
	Reset(Dat);
	For I:=1 To nl Do Begin
		Readln(Dat, E, sigma);
		If (Ei/1000)>E Then Begin {aca se ve cual es la energia anterior a la einc }
			E1:=E;
			s1:=sigma;
			Inc(nm);
		End;
	End;
	Reset(Dat);
	If nm=0 Then begin  {esto es cuando la einc es menor que el primer valor de E}
		Readln(Dat, E1, s1);
		Readln(Dat, E2, s2);
	End
	Else If nm=nl Then begin  {esto es cuando la einc es mayor que el ultimo valor de E}
		For I:=1 To nl-2 Do Readln(Dat);
		Readln(Dat, E1, s1);
		Readln(Dat, E2, s2);
	End
	Else Begin
		For I:=1 To nm Do Readln(Dat);
		Readln(Dat, E2, s2);    {se lee la energia posterior a la einc }
	End;
	E1:=E1*1000;
	E2:=E2*1000;

	Close(Dat);

	SigmaL2:=(((Ei-E1)*(S2-S1))/(E2-E1)+S1)*1e-24;	{de mb a cm^2}

End;

Function SigmaL3(Z:Integer;Ei:Double)  :Double; {datos de PIXE2010: Proton/alpha Ionization (k, l, m shell), Tabulated Cross Section Library.}
	Var
	i,nm,nl				:Integer;
	E1,E2,S1,S2,E,sigma	        :Double;
	Dat				:Text;
        DiSigma                         :String;
Begin
        Case SistOp of
             1: DiSigma := DiUser+'/dat/proton/l/ecpssr/l3-';
             2: DiSigma := DiUser+'\dat\proton\l\ecpssr\l3-';
        End;
	Assign(Dat,DiSigma+FloatToStr(Z)+'.dat');
	Reset(Dat);
	nl:=0; {numero de lineas del archivo}
	nm:=0;
	Repeat		{vemos cuantas lineas tiene el archivo para ese elemento}
		Readln(Dat);
		Inc(nl);
	Until(SeekEOf(Dat));
	nl:=nl-2;
	Reset(Dat);
	For I:=1 To nl Do Begin
		Readln(Dat, E, sigma);
		If (Ei/1000)>E Then Begin {aca se ve cual es la energia anterior a la einc }
			E1:=E;
			s1:=sigma;
			Inc(nm);
		End;
	End;
	Reset(Dat);
	If nm=0 Then begin  {esto es cuando la einc es menor que el primer valor de E}
		Readln(Dat, E1, s1);
		Readln(Dat, E2, s2);
	End
	Else If nm=nl Then begin  {esto es cuando la einc es mayor que el ultimo valor de E}
		For I:=1 To nl-2 Do Readln(Dat);
		Readln(Dat, E1, s1);
		Readln(Dat, E2, s2);
	End
	Else Begin
		For I:=1 To nm Do Readln(Dat);
		Readln(Dat, E2, s2);    {se lee la energia posterior a la einc }
	End;
	E1:=E1*1000;
	E2:=E2*1000;

	Close(Dat);

	SigmaL3:=(((Ei-E1)*(S2-S1))/(E2-E1)+S1)*1e-24;	{de mb a cm^2}

End;

Function SigmaM(Z:Integer;Ei:Double)  :Double; {datos de PIXE2010: Proton/alpha Ionization (k, l, m shell), Tabulated Cross Section Library.}
	Var
	i,nm,nl				:Integer;
	E1,E2,S1,S2,E,sigma	        :Double;
	Dat				:Text;
        DiSigma                         :String;
Begin
        Case SistOp of
             1: DiSigma := DiUser+'/dat/proton/m/ecpssr/m-';
             2: DiSigma := DiUser+'\dat\proton\m\ecpssr\m-';
        End;
	Assign(Dat,DiSigma+FloatToStr(Z)+'.dat');
	Reset(Dat);
	nl:=0; {numero de lineas del archivo}
	nm:=0;
	Repeat		{vemos cuantas lineas tiene el archivo para ese elemento}
		Readln(Dat);
		Inc(nl);
	Until(SeekEOf(Dat));
	nl:=nl-2;
	Reset(Dat);
	For I:=1 To nl Do Begin
		Readln(Dat, E, sigma);
		If (Ei/1000)>E Then Begin {aca se ve cual es la energia anterior a la einc }
			E1:=E;
			s1:=sigma;
			Inc(nm);
		End;
	End;
	Reset(Dat);
	If nm=0 Then begin  {esto es cuando la einc es menor que el primer valor de E}
		Readln(Dat, E1, s1);
		Readln(Dat, E2, s2);
	End
	Else If nm=nl Then begin  {esto es cuando la einc es mayor que el ultimo valor de E}
		For I:=1 To nl-2 Do Readln(Dat);
		Readln(Dat, E1, s1);
		Readln(Dat, E2, s2);
	End
	Else Begin
		For I:=1 To nm Do Readln(Dat);
		Readln(Dat, E2, s2);    {se lee la energia posterior a la einc }
	End;
	E1:=E1*1000;
	E2:=E2*1000;

	Close(Dat);

	SigmaM:=(((Ei-E1)*(S2-S1))/(E2-E1)+S1)*1e-24;	{de mb a cm^2}
       { ShowMessage('SigmaM: '+FloatToStr(SigmaM)); }
End;


Function SigmaK_ion(Z:Integer) :Double; {en cm^2}
	Var
	a0,a1,a2,a3,a4,a5,a6	:Double;
Begin
	If	Z=9 Then Begin 		{Fluor, de 300 a 2000 keV}
		a0:=-6.378278e4;
		a1:= 4.903723e2;
		a2:=-3.303101e-1;		
		a3:=-3.555323e-5;
		a4:= 1.566194e-7;		
		a5:=-7.613083e-11;
		a6:= 1.212132e-14;		
	End
	Else If Z=10 Then Begin	{Neon, de 100 a 1200 keV}
		a0:= 6.848821e3;
		a1:=-1.610356e2;
		a2:= 1.300972e0;	
		a3:=-2.677139e-3;
		a4:= 2.779273e-6;		
		a5:=-1.485295e-09;
		a6:= 3.232778e-13;		
	End	
	Else If	Z=12 Then Begin	{Magnesio, de 100 a 1700 keV}
		a0:= 1.066415e3;
		a1:=-2.484061e1;
		a2:= 1.767489e-1;		
		a3:=-2.259626e-4;
		a4:= 1.406597e-7;		
		a5:=-4.452442e-11;
		a6:= 5.704357e-15;		
	End			
	Else If Z=13 Then Begin	{Aluminio, de 100 a 3000 keV}
		a0:= 6.673815e2;
		a1:=-1.340031e1;
		a2:= 8.013209e-2;		
		a3:=-7.996898e-5;
		a4:= 3.751489e-8;		
		a5:=-8.703737e-12;
		a6:= 8.015777e-16;	
	End
	Else If	Z=14 Then Begin	{Silicio, de 25 a 4000 keV}
		a0:= 6.921673e1;
		a1:=-2.948013e0;
		a2:= 2.905764e-2;		
		a3:=-2.294673e-5;
		a4:= 8.185490e-9;		
		a5:=-1.424149e-12;
		a6:= 9.770609e-17;		
	End
	Else If	Z=15 Then Begin	{Fosforo, de 600 a 2600 keV}
		a0:= 2.308992e3;
		a1:=-1.577773e1;
		a2:= 4.487383e-2;		
		a3:=-3.931785e-5;
		a4:= 1.773023e-8;		
		a5:=-4.119134e-12;
		a6:= 3.881263e-16;		
	End
	Else If Z=16 Then Begin	{Azufre, de 600 a 5000 keV}
		a0:=-1.267672e3;
		a1:= 2.254964e0;
		a2:= 5.427487e-3;		
		a3:=-3.493528e-6;
		a4:= 9.275580e-10;	
		a5:=-1.179757e-13;
		a6:= 5.895898e-18;		
	End	
	Else If Z=17 Then Begin	{Cloro , de 800 a 5000 keV}
		a0:=-1.258650e3;
		a1:= 2.170627e0;
		a2:= 2.656058e-3;		
		a3:=-1.581386e-6;
		a4:= 3.676792e-10;	
		a5:=-3.974883e-14;
		a6:= 1.620375e-18;		
	End		
	Else If	Z=20 Then Begin	{Calsio, de 500 a 2000 keV}
		a0:= 1.154384e2;
		a1:=-6.917885e-1;
		a2:= 1.333803e-3;		
		a3:= 1.889016e-7;
		a4:=-5.930152e-10;		
		a5:= 2.591265e-13;
		a6:=-3.875577e-17;		
	End
	Else Writeln('No esta modelada la seccion eficaz para este elemento');
	
	SigmaK_ion:=(a0+a1*Einc+a2*Einc*Einc+a3*Einc*Einc*Einc+a4*Einc*Einc*Einc*Einc+a5*Einc*Einc*Einc*Einc*Einc+
				a6*Einc*Einc*Einc*Einc*Einc*Einc)*1e-24; {para pasar de barn a cm^2}
End;	


Procedure Refina(); 
	LABEL 99;
	Var
	NPR,i,j,u,l,NV,m	:Integer;
	Sal          		:Text;
        En           		:Double;
        PPP			:pMN;
        PP		:DoubleArrayNP;
Begin
        {TablaEfi();      }
        If NFilm=0 Then TablaSP();
       {Guardamos la seccion eficaz de las lineas}
        If NFilm<>0 Then Begin
            For u:=1 To Ndet+Nfdet Do For j:=1 To Nlin[u] Do Begin
                Case Borde[u,j][1] Of
                      'K': S[u,j]:=SigmaK(Zdet[u],Einc);
                      'L':Begin
                                 Case Borde[U,J][2] Of
                                 '3': S[u,j]:=SigmaL3(Zdet[u],Einc)+CK[Zdet[U],3]*SigmaL2(Zdet[u],Einc)+(CK[Zdet[U],2]+CK[Zdet[U],1]*CK[Zdet[U],3])*SigmaL1(Zdet[u],Einc);
                                 '2': S[u,j]:=SigmaL2(Zdet[u],Einc)+CK[Zdet[U],1]*SigmaL1(Zdet[u],Einc);
                                 '1': S[u,j]:=SigmaL1(Zdet[u],Einc);
                                 End
                      End;
                      'M': S[u,j]:=SigmaM(Zdet[u],Einc);
                End;
            End;
        End
        Else Begin
             Nint:=Round(Einc/10.0);
             For i:=1 To Nint Do Begin
                 Eint[i]:=i*Einc/Nint;
                 For u:=1 To Ndet Do For j:=1 To Nlin[u] Do Begin
                     Case Borde[u,j][1] Of
                          'K': CS[i,u,j]:=SigmaK(Zdet[u],Eint[i]);
                          'L':Begin
                               Case Borde[U,J][2] Of
                                    '3': CS[i,u,j]:=SigmaL3(Zdet[u],Eint[i])+CK[Zdet[U],3]*SigmaL2(Zdet[u],Eint[i])+(CK[Zdet[U],2]+CK[Zdet[U],1]*CK[Zdet[U],3])*SigmaL1(Zdet[u],Eint[i]);
                                    '2': CS[i,u,j]:=SigmaL2(Zdet[u],Eint[i])+CK[Zdet[U],1]*SigmaL1(Zdet[u],Eint[i]);
                                    '1': CS[i,u,j]:=SigmaL1(Zdet[u],Eint[i]);
                               End
                          End;
                          'M': Begin
                               CS[i,u,j]:=SigmaM(Zdet[u],Einc);
                          End;
                     End;
                 End;
             End;
        End;
        {Creando el arreglo P_ref de parametros que se van a refinar y el arreglo de lam_ref correspondiente a los parametros que se refinan}
	For i:=1 To NPG+4*NL+NG+Nel+4*NLF+Nfilm Do Begin
		Sumcod[i]:=Sumcod[i-1]+Cod[i];
		NPR:=Sumcod[i];{numero de parametros a refinar}
		If Cod[i]=1 Then Begin
			P_Ref[Sumcod[i]]:=P0[i];
			Lam_Ref[Sumcod[i]]:=Lam[i];
		End;	  	
	End;
	{Para cuando no se ajusta ningun parametro}
	If NPR=0 Then GOTO 99;

	{Esto nos sirve para la condicion de que la suma de los factores de escala del mismo grupo sea 1}
	For i:=0 To 60 Do Begin
		GrupoV[i]:=0; {contiene el Numero de Grupo vinculado con la linea I (solo las que se refinan)}
		Sref[i]:=0;	  {suma de las fracciones refinadas del grupo I} 
		SnoRef[i]:=0  {suma de las fracciones no refinadas del grupo I}
	End;
	NV:=0;{numero de vinculo, ej: 1->kk, 2->l3, etc}
	l:=0;
	For i:=1 To Ndet Do Begin
		Inc(NV);
		For j:=1 To Nlin[i] do Begin
			Inc(l);
			If j>1 Then If (Borde[i,j]<>Borde[i,j-1]) Then Inc(NV);
			If Cod[NPG+NL+l]=1 Then Begin
				GrupoV[Sumcod[NPG+NL+l]]:=NV;
				Sref[NV]:=Sref[NV]+P_ref[Sumcod[NPG+NL+l]]
			End
			Else SNoref[NV]:=SNoref[NV]+P0[NPG+NL+l]
		End
	End;
	Inc(NV);  {Las concentraciones las agregamos como un vinculo mas para normalizarlas junto con las fracciones}
	For i:=1 To Nel Do Begin
		If Cod[NPG+NG+4*NL+i]=1 Then Begin
			GrupoV[Sumcod[NPG+NG+4*NL+i]]:=NV;
			Sref[NV]:=Sref[NV]+P_ref[Sumcod[NPG+NG+4*NL+i]]
		End
		Else SNoref[NV]:=SNoref[NV]+P0[NPG+NG+4*NL+i]
	End;
	{Ahora hacemos lo mismo para el film}
	l:=0;	
	For i:=Ndet+1 To Ndet+Nfdet Do Begin
		Inc(NV);
		For j:=1 To Nlin[i] do Begin
			Inc(l);
			If j>1 Then If (Borde[i,j]<>Borde[i,j-1]) Then Inc(NV);
			If Cod[NPG+NG+4*NL+Nel+NLF+l]=1 Then Begin
				GrupoV[Sumcod[NPG+NG+4*NL+Nel+NLF+l]]:=NV;
				Sref[NV]:=Sref[NV]+P_ref[Sumcod[NPG+NG+4*NL+Nel+NLF+l]]
			End
			Else SNoref[NV]:=SNoref[NV]+P0[NPG+NG+4*NL+Nel+NLF+l]
		End
	End;
	Inc(NV);  
	For i:=1 To Nfilm Do Begin
		If Cod[NPG+NG+4*NL+Nel+4*NLF+i]=1 Then Begin
			GrupoV[Sumcod[NPG+NG+4*NL+Nel+4*NLF+i]]:=NV;
			Sref[NV]:=Sref[NV]+P_ref[Sumcod[NPG+NG+4*NL+Nel+4*NLF+i]]
		End
		Else SNoref[NV]:=SNoref[NV]+P0[NPG+NG+4*NL+Nel+4*NLF+i]
	End;
	{Normalizacion de las fracciones de linea}
	For i:=1 To NPR Do If GrupoV[i]<>0 Then P_ref[i]:=P_ref[i]*(1.0-Snoref[GrupoV[i]])/Sref[GrupoV[i]];

	{Creando la matriz que contiene los NPR+1 conjuntos de NPR parametros}
	For j:=1 To NPR Do PPP[1,j]:=P_Ref[j];  

	For j:=2 To NPR+1 Do Begin
        For i:=1 to NPR do Sref[GrupoV[i]]:=0.0;
		For i:=1 To NPR Do Begin
			If i<>j-1 Then PP[i]:=P_Ref[i]
			Else If (cod[2]=1) and (I=sumcod[2]) Then PP[i]:=P_Ref[i]+Lam_Ref[i]
            Else PP[i]:=P_Ref[i]+Lam_Ref[i]*P_Ref[i];
            If GrupoV[i]<>0 Then Sref[GrupoV[i]]:=Sref[GrupoV[i]]+PP[i] {para luego normalizar las fr}
		End;	
		For i:=1 To NPR Do Begin
			If GrupoV[i]<>0 Then Begin
				PP[i]:=PP[i]*(1.0-Snoref[GrupoV[i]])/Sref[GrupoV[i]]; {normalizacio de las fr}
			End;
			PPP[j,i]:=PP[i];
		End;	
	End;

	{Refinando los parametros con Amoeba}
	P_final:=Amoeba(PPP,NPR);

	{Escribimos el P0 con los nuevos parametros}	
	For i:=1 To NPG+NG+4*NL+Nel+4*NLF+Nfilm Do Begin
		If Sumcod[i]>0 Then P0[i]:=(1.0-Cod[i])*P0[i]+Cod[i]*P_final[Sumcod[i]]
                Else P0[i]:=P0[i]
	End;
	For i:=1 To Nel Do C[i]:=P0[NPG+NG+4*NL+i];
	l:=0;
	For i:=Nel+1 To Nel+Nfilm Do Begin
		Inc(l);
		C[i]:=P0[NPG+NG+4*NL+Nel+4*NLF+l];
	End;	
	l:=0;
	For i:=1 To Ndet Do For j:=1 To Nlin[i] Do Begin
		Inc(l);
		E[i,j]:=P0[NPG+l];
		Fr[i,j]:=P0[NPG+NL+l];
                Tail[i,j]:=P0[NPG+2*NL+l];
                Beta[i,j]:=P0[NPG+3*NL+l];
	End;
	l:=0;
	For i:=Ndet+1 To Ndet+Nfdet Do For j:=1 To Nlin[i] Do Begin
		Inc(l);
		E[i,j]:=P0[NPG+NG+4*NL+Nel+l];
		Fr[i,j]:=P0[NPG+NG+4*NL+Nel+NLF+l];
                Tail[i,j]:=P0[NPG+NG+4*NL+Nel+2*NLF+l];
                Beta[i,j]:=P0[NPG+NG+4*NL+Nel+3*NLF+l];
	End;	
	For i:=1 To NG Do Fesc[i]:=P0[NPG+NG+4*NL+i];


       { Errores;
        For I:=1 To NPG+NG+4*NL+Nel+4*NLF+Nfilm Do If Cod[I]=1 Then PE[I]:=PE_ref[SumCod[I]]; }
	99:

        If(NPR=0) Then Begin
             Chi2:=Func(P_final,NPR);
             Form1.label2.Caption:=FloatToStr(Iter);
             Form1.label4.Caption:=FloatToStr(Chi2);
             Form1.Refresh;
         End;
End;

Procedure EscribirSal;
   Var
       Sal  :Text;
       i,j,m,u    :Integer;
       En     :Double;
Begin
{Escribiendo el archivo .sal para graficar}
Assign(Sal,DiSal+'.sal');
Rewrite(Sal);
For i:=Canali to Canalf Do Begin
En:=(P0[1]*I+P0[2]);
Write(Sal,En:7:4,'  ',Ymed[i],'   ',Ycal[i]:10:4,'   ',Yfondo[i]:10:4);
For u:=1 To Ndet+Nfdet Do For j:=1 To Nlin[u] Do Write(Sal,' ',Ylinea[i,u,j]:10:4);
For u:=1 To Ndet+Nfdet Do For j:=1 To Nlin[u] Do Write(Sal,' ',Ycola[i,u,j]:10:4);
For m:=1 To Npicosuma Do Write(Sal,' ',(ISuma[m]*X[11])*Perfil_RAS(En,Esuma[m],Sigmasuma2[m]):10:4);
Writeln(Sal);
End;
Close(sal);

End;


Function Amoeba(PPP :pMN; NPR :Integer) :DoubleArrayNP;
	Var
	alpha,beta,gamma,sum    :Double;
	YY						:DoubleArrayMP;
	i,j,h,l,mayor,K	:Integer;
	Pc,P1,P2				:DoubleArrayNP;
	Y1,Y2,Dif			:Double;
Begin
	alpha:=1.0;{factor de reflexion}
	beta:=0.5; {factor de contraccion}
	gamma:=2.0;{factor de expansion}
	l:=0;
        YY[l]:=0.0;
	iter:=0;
	Dif:=100;
	While Dif>Titer Do Begin {sigue hasta que Dif sea menor que Titer por mas de 19 iteraciones seguidas}
		Inc(iter);
        For i:=1 to NPR do Sref[GrupoV[i]]:=0.0;
		{YY es el arreglo de la funcion func evaluada en los distintos grupos de parametros}
		For i:=1 To MP Do YY[i]:=0.0; 
		For i:=1 To NPR+1 Do YY[i]:=Func(PPP[i],NPR);
		{Determinamos cual grupo de parametros da un Y mayor(h) y cual da un Y menor(l)}
		h:=1;
		l:=1;
		For i:=2 To NPR+1 Do Begin
			If YY[i]>YY[h] Then h:=i;
			If YY[i]<YY[l] Then l:=i;
		End;
		Dif:=2.0*abs(YY[h]-YY[l])/(abs(YY[h])+abs(YY[l]));

                Chi2:=Func(PPP[l],NPR);
                Form1.label2.Caption:=FloatToStr(Iter);
                Form1.label4.Caption:=FloatToStr(Chi2);
                Form1.Refresh;

		{Calculamos el centroide Pc de todos los grupos de parametros}
		For i:=1 To NPR Do Begin
			sum:=0.0;
			For j:=1 To NPR+1 Do sum:=sum+PPP[j,i];
			Pc[i]:=sum/(NPR+1);
		End;
		
		{comienza el amoeba}
		For i:=1 To NPR Do P1[i]:=(1+alpha)*Pc[i]-alpha*PPP[h,i]; 	{Reflexion}
		Y1:=Func(P1,NPR);
		If Y1<YY[l] Then Begin
			For i:=1 To NPR Do P2[i]:=(1+gamma)*P1[i]-gamma*Pc[i];   {Expansion}
			Y2:=Func(P2,NPR);
			If Y2<YY[l]	Then PPP[h]:=P2
			Else PPP[h]:=P1;
		End	
		Else Begin
			mayor:=0;
			For i:=2 To NPR+1 Do Begin {Vemos si y1 es mayor que todos los yy[i] con i distinto de h}
				If i<>h Then If Y1>YY[i] Then Inc(mayor);
			End;
			If mayor=NPR-1 Then Begin
				If Y1<YY[h] Then PPP[h]:=P1;
				For i:=1 To NPR Do P2[i]:=beta*PPP[h,i]+(1-beta)*Pc[i];  {Contraccion}
                                Y2:=Func(P2,NPR);
				If Y2>YY[h] Then Begin
					For i:=1 To NPR+1 Do Begin 
						For j:=1 To NPR Do Begin 
							PPP[i,j]:=(PPP[i,j]+PPP[l,j])/2.0;
							If GrupoV[j]<>0 Then Sref[GrupoV[j]]:=Sref[GrupoV[j]]+PPP[i,j]; {normalizacion de fr}
						End;	
						For j:=1 To NPR Do If GrupoV[j]<>0 Then Begin
							PPP[i,j]:=PPP[i,j]*(1.0-SNoref[GrupoV[j]])/Sref[GrupoV[j]];	
						End;		
						For k:=1 To NPR Do Sref[GrupoV[k]]:=0.0;
					End;
				End	
				Else PPP[h]:=P2; 
			End  
			Else PPP[h]:=P1;
		End;
		{normalizacion de fr}
		For j:=1 To NPR Do Sref[GrupoV[j]]:=0.0;
		For j:=1 To NPR Do If GrupoV[j]<>0 Then Sref[GrupoV[j]]:=Sref[GrupoV[j]]+PPP[h,j];
		For j:=1 To NPR Do If GrupoV[j]<>0 Then Begin
			PPP[h,j]:=PPP[h,j]*(1.0-SNoref[GrupoV[j]])/Sref[GrupoV[j]];
		End;	
	End;
        {ShowMessage('Final Amoeba');}
	Amoeba:=PPP[l];
End;


Procedure Picosuma;
	Var
	Ipico			:Array[1..30,1..30] Of Double;
	L,U,J,UU,JJ 		:Integer;
	Isum,Sigmasum2,Esum  	:Double;
Begin

   L:=0;
   If NFilm=0 Then Begin
	For U:=1 To Ndet Do For J:=1 To Nlin[U] Do Begin
		Suma[U,J]:=0.0;
                Ipico[U,J]:=Intensidad_mg2(U,J,E[U,J]);
                For UU:=1 To Ndet Do For JJ:=1 To Nlin[UU] Do Begin
                      Ipico[UU,JJ]:=Intensidad_mg2(UU,JJ,E[UU,JJ]);
                      ISum:=Ipico[U,J]*Ipico[UU,JJ];					 {Intensidad del pico suma}
		      Esum:=E[U,J]+E[UU,JJ]; 						 {Energia del pico suma}
		      Sigmasum2:=X[9]*X[9]+0.00376*X[10]*Esum;	 			 {Ancho del pico suma}
		      If (Esum<Einc) Then
		      If ((Isum*X[11])>0.1*Sqrt(Sqrt(SigmaSum2)*Efi(E[U,J]+E[UU,JJ])*Fondo(E[U,J]+E[UU,JJ],X[6],X[7],X[8]))) Then Begin
				Inc(L);				{primero verificamos que cumpla las condiciones y luego guardamos la informacion que necesitamos}
				Isuma[L]:=Isum;
				Sigmasuma2[L]:=Sigmasum2;
				Esuma[L]:=Esum;
				If (U<>UU) Or (J<>JJ) Then Begin
					Suma[U,J]:=Suma[U,J]+Ipico[UU,JJ]/2.0;
					Isuma[L]:=Isum/2.0 {se divide por 2 por que para picos distintos se crean dos picos sumas con igual posicion e intensidad}
				End
				Else Suma[U,J]:=Suma[U,J]+Ipico[UU,JJ];
			End;			
		End;				
	End;
   end
   Else Begin
	For U:=1 To Ndet+Nfdet Do For J:=1 To Nlin[U] Do Begin
		Suma[U,J]:=0.0;
		If U<=Ndet Then Ipico[U,J]:=Intensidad(U,J,E[U,J],S[U,J])	 {Intensidad del primer pico(tenemos que ver si es del}
		Else If U>Ndet Then Ipico[U,J]:=Intensidad_film(U,J,E[U,J],S[U,J]); {sustrato o del film)}
                For UU:=1 To Ndet+Nfdet Do For JJ:=1 To Nlin[UU] Do Begin
		      If UU<=Ndet Then Ipico[UU,JJ]:=Intensidad(UU,JJ,E[UU,JJ],S[UU,JJ]) {Intensidad del primer pico(tenemos que ver si es del}
		      Else If UU>Ndet Then Ipico[UU,JJ]:=Intensidad_film(UU,JJ,E[UU,JJ],S[UU,JJ]); {sustrato o del film)}
                      ISum:=Ipico[U,J]*Ipico[UU,JJ];					 {Intensidad del pico suma}
		      Esum:=E[U,J]+E[UU,JJ]; 						 {Energia del pico suma}
		      Sigmasum2:=X[9]*X[9]+0.00376*X[10]*Esum;	 			 {Ancho del pico suma}
		      If (Esum<Einc) Then
		      If ((Isum*X[11])>0.1*Sqrt(Sqrt(SigmaSum2)*Efi(E[U,J]+E[UU,JJ])*Fondo(E[U,J]+E[UU,JJ],X[6],X[7],X[8]))) Then Begin
				Inc(L);				{primero verificamos que cumpla las condiciones y luego guardamos la informacion que necesitamos}
				Isuma[L]:=Isum;
				Sigmasuma2[L]:=Sigmasum2;
				Esuma[L]:=Esum;
				If (U<>UU) Or (J<>JJ) Then Begin
					Suma[U,J]:=Suma[U,J]+Ipico[UU,JJ]/2.0;
					Isuma[L]:=Isum/2.0 {se divide por 2 por que para picos distintos se crean dos picos sumas con igual posicion e intensidad}
				End
				Else Suma[U,J]:=Suma[U,J]+Ipico[UU,JJ];
			End;
		End;
	End
   end;

   Npicosuma:=L;
End;

Function Resc(U,J:Integer):Double; {Pico de escape}
Var
 Tita:Double;
Begin
 Tita:=0.0;
 If (E[U,J]>1.84) Then
     Resc:=0.0202/(1.0+((0.01517*Cos(Tita)-0.000803)*E[U,J]+0.0455*Cos(Tita)+0.01238)*Sqr(E[U,J]))
 Else Resc:=0.0      {Scott, Love y Reed 2nd ed. p\A0g. 78}
End;


End. {fin del archivo}
