program SimFitEAv;
     var i, j, k, nc, pt, datawindowID, maxsamples, r,
         drawwindowID, prleft, prtop, prright, prbottom, ltrunc, rectH, dlogFWindowID, ihalf :integer;
 
         bg, os, l, lowlim, hilim, EA, F0, b, bglow: array[1..16] of real;
         dF, dFsum, Fs: array[1..50] of real;
                            
 		 DWname, ns, eas, s:string;
 
 		 dFref, w, EF, dFsummax, maxfit, cscale, min, cbottom, ctop, c, maxspan, fl :real;
 
begin
  maxsamples := 16;
  
  writeln('new run');
         
  SetBoxTitle('New Data Window');
  DWname :='';
  Input('Name of data window?', DWname);    
  NewDataWindow(name DWname + ' Data', nrRows 100, nrCols 2);
  datawindowID := GetCurrentWindow(dataType);
 
 
  Alert('Copy data (max. 16 curves) from Excel, then click OK');
  Paste;
 
  nc:=nrcols/2;
  
  i:=0;   
  repeat  
    i := i+1;
  until DataOK(i+1, 1) = false;
  pt := i; {points per curve}
     
  InsertColumns(at nrcols+1, count 4);
 
    
      
  for i := 1 to nc do
    begin
      SetColumnProperties(col i*2-1, name 'Cycle');
      NumberToString(i,ns,true,0);
 	  SetColumnProperties(col i*2, name 'Fl Curve ' + ns);
    end;
    
  SetDataWindowProperties(window datawindowID, nrrows 100 * nc);
   
  SetColumnProperties(col nrcols-1, name 'Cycle');
  SetColumnProperties(col nrcols, name 'Selected Points');
  
  SetColumnProperties(col nrcols-3, name 'sim cycle');
  SetColumnProperties(col nrcols-2, name 'dF');
 
  for j:= 1 to nc do {curve}
    for i:= 2 to pt do {cycle}
      begin
        data[(j-1)*50 + i, nrcols-3] := i;
        data[(j-1)*50 + i, nrcols-2] := data[i, j*2] - data[i-1, j*2]; {prepare dF data}
      end;  
  
  
  
  
  j := 1; {Kurve #}
  
  repeat
    
    {find highest dF 5 point peak;
      5 point peak with lower rel. err. than 4 point or 3 point peak}
    
    for i := 2 to pt do dF[i] := data[(j-1) * 50 + i, nrcols - 2];
    for i := 8 to pt do dFsum[i] := dF[i - 4] + dF[i - 3] + dF[i - 2] + dF[i - 1] + dF[i]; {alle mšglichen 5 point peaks summieren} 
    dFsummax := dFsum[8];
    for i := 9 to pt do
      if dFsum[i] > dFsummax then
          begin
            dFsummax := dFsum[i];
            l[j] := i; {der letzte Datenpunkt des hšchsten 5 point peaks ist das Limit; spŠtere Zyklen gehen nicht in die Auswertung ein}
          end;
       
    {jetzt rŸckwŠrts gehen und die vorlaufende Grundlinie aus 9 Punkten festlegen;
     Kriterium: 3 Punkte liegen "hšher" als Mittel von 3 Testpunkten}
    
    i := l[j] - 4;
    repeat
      i := i - 1;
      dFref := (dF[i] + dF[i+1] + dF[i+2])/3; {reference level}
      w := 0;
      for k := 1 to 8 do if (dF[i-k] > dFref) then w := w + 1;
    until (w > 2) or (i = 9); 
    hilim[j] := i;
    
          
    {Einzelanpassung machen}
    
    SelectFunction('EAvPeak');
    
    w := 0.5 * (data[pt, j*2] - data[1, j*2]); {nŠherungsweise 0.5 * F curve amplitude}
    fl := data[(j-1) * 50 + l[j], nrcols - 2]; {F des letzten brauchbaren Punktes} 
    fl := fl / (1 + 0.9/exp(fl^2/w^2))^l[j];    
       
    SetParameterProperties(param 1, value 0.9); {amp}
    SetParameterProperties(param 2, value fl, min 0); {F0}
    SetParameterProperties(param 3, value w^2, min 0); {b}
    SetParameterProperties(param 4, value 0); {background slope}
        
    selectwindow(datawindowID);
    SelectCell(fromRow (j-1)*50 + l[j] + 1, toRow j*50, fromCol nrcols-3, toCol nrcols-2); {clear unsuitable data points}
    Clear;
    SelectCell(fromRow (j-1)*50 + hilim[j] - 8, toRow (j-1)*50 + l[j], fromCol nrcols-3, toCol nrcols-2);
    
    Fit(algorithm levenberg, xColumn nrcols-3, yColumn nrcols-2, xErrKind zeroError, yErrKind unknownError, yError 1, selRowsOnly true,
    printResults false, onlyActiveParams true);
    Domenu('Calc:Get Fitted Params');
        
    EA[j] := GetResult(fittedParameter, 1); {EA}
    F0[j] := GetResult(fittedParameter, 2); {F0}
    b[j] := GetResult(fittedParameter, 3); {peak amplitude}
    bg[j] := GetResult(fittedParameter, 4); {background slope}
    
    
    
    {jetzt mit F-Daten den Background bestimmen}
        
    if hilim[j] < 9 then hilim[j] := 9; {Absicherung; falls der Peak schlecht definiert ist und die Nulllinie dann hšher liegt als alle Geradenpunkte}
    lowlim[j] := hilim[j] - 8;
      
    SelectCell(fromRow lowlim[j], toRow hilim[j], fromCol j*2-1, toCol j*2);
    SelectFunction('Polynomial');   
    SetParameterProperties(param 1, value 1); {Gerade}
    SetParameterProperties(param 2, value 2); {offset}
    SetParameterProperties(param 3, value 0); {slope}
    Fit(algorithm linear, xColumn j*2-1, yColumn j*2, xErrKind zeroError, yErrKind unknownError, selRowsOnly true,
    printResults false, onlyActiveParams true); {-99: use row index} 
    Domenu('Calc:Get Fitted Params');
    bg[j] := GetResult(fittedParameter, 3); {background slope}
    os[j] := GetResult(fittedParameter, 2); {background offset}
    
    
    
    {background subtrahieren, Daten tabellieren}

    setcurrentwindow(datawindowID);    
    selectwindow(datawindowID);
         
    SelectCell(fromRow lowlim[j], toRow l[j], fromCol j*2-1, toCol j*2); {transfer selected points}
    copy; 
    selectcell(row lowlim[j] + (j-1) * 50, col nrcols-1);
    paste;
    
    for i:= lowlim[j] to l[j]  do data[i + (j-1) * 50 , nrcols] := data[i + (j-1) * 50, nrcols] - os[j] - bg[j]*i; {subtract background}
        
        
    {jetzt einzeln fitten; kein background mehr}    
           
    SelectFunction('EAv'); 
    SetParameterProperties(param 1, value EA[j], mode paramActive); {EA}
    SetParameterProperties(param 2, value F0[j], min 0, mode paramActive); {f0}
    SetParameterProperties(param 3, value b[j], min 0, mode paramActive); {b}
    SetParameterProperties(param 4, value 0, mode paramInactive); {background slope}
    SetParameterProperties(param 5, value 0, mode paramInactive); {background offset}
            
    SelectCell(fromRow lowlim[j] + (j-1) * 50, toRow l[j] + (j-1) * 50, fromCol nrcols-1, toCol nrcols);
    
    Fit(algorithm robust, xColumn nrcols-1, yColumn nrcols, xErrKind zeroError, yErrKind unknownError, selRowsOnly true,
    printResults false, onlyActiveParams true); {-99: use row index} 
    Domenu('Calc:Get Fitted Params');
    
    EA[j] := GetResult(fittedParameter, 1); {EA}
    F0[j] := GetResult(fittedParameter, 2); {f0}
    b[j] := GetResult(fittedParameter, 3); {b}
      
    writeln('EA: ', 1 + EA[j]);
    
    for i:= lowlim[j] + (j-1) * 50 to l[j] + (j-1) * 50 do data[i, nrcols-1] := data[i, nrcols-1] + (j-1) * 50; {Zeilennummern fŸr simultanes fitten umstellen}
       
       
    j := j + 1;
  until j > nc;       
  
  
  
        
        
  {#########################################################################}
  {################     simultanes fitting         #########################}
  {#########################################################################}
  
  
  EF := 0;
  for i := 1 to nc do EF := EF + EA[i];
  EF := EF/nc; {average EF for simult. fitting}
    
  SetCurrentWindow(datawindowID);
  SelectWindow(datawindowID);

  
  {SetBoxTitle('Simultaneous Fitting of Real-Time PCR Data');
  fixEF := 0;
  Input('$WSelect data window',datawindowID, 'Fixed EF: ', fixEF); }

  
  SetDefaultCols (nrcols-1, nrcols, 0, 0);
  
  SelectFunction('M16EAv');
  for i:= 1 to nc do SetParameterProperties(param i, value b[i], min 0, mode paramActive); {startwert b aus individ. Kurvenfit}
  for i:= 1 to nc do
    begin
      c := l[i]-2;
      c := F0[i]*(1+EA[i])^c/(1+EF)^c; {F0 an average EF anpassen, damit fitting besser funktioniert}   
      SetParameterProperties(param 16 + i, value c, min 0, mode paramActive); {startwert F0}
    end;
         
  if nc < 16 then
    for i:= nc+1 to 16 do
      begin
        SetParameterProperties(param i, value 0, mode paramUnused);
        SetParameterProperties(param 16+i, value 0, mode paramUnused);
      end;
 
  SetParameterProperties(param 33, value EF, mode paramActive); {startwert EA}
     
  Fit(algorithm levenberg, xColumn nrcols-1, yColumn nrcols, xErrKind zeroError, yErrKind unknownError, selRowsOnly false,
    printResults true, onlyActiveParams true);
  Domenu('Calc:Get Fitted Params');
   
  Fit(algorithm levenberg, xColumn nrcols-1, yColumn nrcols, xErrKind zeroError, yErrKind unknownError, selRowsOnly false,
    printResults true, onlyActiveParams true);
  Domenu('Calc:Get Fitted Params'); {Doppel-Fit: hŠsslich, aber in EinzelfŠllen hilfreich}

  EF := GetResult(fittedParameter, maxsamples * 2 + 1);
  
 
 
       
        
 {#########################################################################}
 {#############             Resultate ausgeben               ##############}
 {#########################################################################}


  s := GetWindowProperty(datawindowID, name);
  Delete(s,length(s)-4,7); {"Data" abziehen}
  NewDrawingWindow(name  s + ' Draw');
    
  drawwindowID := GetCurrentWindow(drawingType);
  DoMenu('File:Page Setup...');
      
  SetTextStyle('Helvetica Regular',18,bold);
  MoveTo(prleft + 90, prtop + 10);
  DrawText(s,0,false); 
  
  {max. Spanne suchen}
  
  maxspan := 0;
  for j:= 1 to nc do {jede Kurve einzeln auswerten}
    begin
      
      maxfit := 0; {max. raw Fluoreszenz eines fit-Punktes suchen}
      for i:= 1 to l[j] do if DataOK((j-1)*50 + i, nrcols) = true then
        begin
          w := data[i, j*2];
          if w > maxfit then maxfit := w;
        end;    
   
      min := 100; {min. raw Fluoreszenz suchen}
      for i:= 1 to l[j] do if DataOK(i, j*2) = true then
        begin
          w := data[i, j*2];
          if w < min then min := w;
        end; 
      bglow[j] := min;
      
      if (maxfit - min) > maxspan then maxspan := (maxfit - min);     
    end;
 
  cscale := maxspan / 0.75; {= 100%; d.h. die hšchste Spanne macht 75% der y-Spanne aus}
  
  prleft := GetWindowProperty(drawwindowID, pageRectLeft);
  prtop := GetWindowProperty(drawwindowID, pageRectTop);
  prright := GetWindowProperty(drawwindowID, pageRectRight);
  prbottom := GetWindowProperty(drawwindowID, pageRectBottom);
  rectH := prbottom - (30 * maxsamples + 10) - prtop; 
      
  for i := 1 to nc do
    begin
            
      b[i] := GetResult(fittedParameter, i); {b-Wert}
      F0[i] := GetResult(fittedParameter, maxsamples + i); {F0 value}
                                        
      k:=1;
      for j := 1 to nc do if GetResult(fittedParameter, maxsamples + j) > F0[i] then k := k + 1; {k: pos. in F0-Abfolge}   
            
      cbottom := bglow[i] - cscale * 0.05;
      ctop := cbottom + cscale;
                        
      w := prtop + 10 + (k-1)*30;         
      
      SetNewGraphRect(prleft+90, w, prright-150, w + rectH);
      
      PlotData(plotType scatterPlot, xColumn i*2-1, yColumn i*2,  window datawindowID, newGraph true, yScaling linScaling,
        connected false, pointType 7, pointSize 4, pointThickness 0.5, pointsRed 0x0000, pointsGreen 0x0000, pointsBlue 0x0000,
        legendEntry false, xlabel '', ylabel '', title '', fillPattern 2);  
    
      
      ClearTicks(xaxis);
      SetAxisProperties(axis 'X1', first 0, last pt, drawLine false, drawLabels false);
      ClearTicks(yaxis);
      SetAxisProperties(axis 'Y1', first cbottom, last ctop, drawLine false, drawLabels false);
      if (ObjectExists(GetAxisObject(axis 'X2'))) <> 0 then DeleteAxis('X2');
      if (ObjectExists(GetAxisObject(axis 'Y2'))) <> 0 then DeleteAxis('Y2');
      
      if k = nc then
        begin
          SetAxisAttributes(xAxis, MinusSideticks);
          SetAxisProperties(axis 'X1', position cbottom, drawLine true, drawLabels true, labelssize 12, labelsstyle plain);
                
          MakeTicks(xaxis,10,10,0);
          
          SetTextStyle('Helvetica Regular',14,plain);
          MoveTo(prleft + 90 + (prright-prleft-240)/2, w + rectH + 32);
          DrawText('Cycle number',0,true);          
        
          SetAxisAttributes(yAxis, MinusSideticks);
          SetAxisProperties(axis 'Y1', position 0, drawLine true, drawLabels true, labelssize 12, labelsstyle plain);
          
          MoveTo(prleft + 50, prtop+10+(k-1)*30 + c + rectH/2);
          DrawText('Fluorescence (bottom curve)',90,true);      
          
        end;
      
      {mark used amplification points}
      
      SelectCell(fromRow hilim[i] + 1, toRow l[i], fromCol i*2-1, toCol i*2);
                       
      PlotData(plotType scatterPlot, xColumn i*2-1, yColumn i*2,  selRowsOnly true, window datawindowID, 
        connected false, pointType 6, pointSize 3, pointThickness 0.5, pointsRed 0x7777, pointsGreen 0x7777, pointsBlue 0x7777,
        legendEntry false, fillPattern 2);
            
      {mark used line points}
            
      SelectCell(fromRow lowlim[i], toRow hilim[i], fromCol i*2-1, toCol i*2);
                      
      PlotData(plotType scatterPlot, xColumn i*2-1, yColumn i*2,  selRowsOnly true, window datawindowID, 
        connected false, pointType 6, pointSize 3, pointThickness 0.5, pointsRed 0x0000, pointsGreen 0x0000, pointsBlue 0x0000,
        legendEntry false, fillPattern 2);
                     
      {display fit result numbers}               
      
      MoveTo(prright-140, 20 + (k-1) * 30);
      NumberToString(1 + EA[i],eas,false,-3);
      NumberToString(F0[i],s,false,3);
      
      SetTextStyle('Helvetica Regular',12,plain);
      DrawText(GetColumnProperty(i*2, name) + '  ea: ' + eas + Chr(13) + 'F0: ' + s , 0, false);
      
      
      {plot fitted function}
       
      SelectFunction('EAv');
      SetFunctionParam('', 1, EF);    {EA}
      SetFunctionParam('', 2, F0[i]); {f0}
      SetFunctionParam('', 3, b[i]);  {b}
      SetFunctionParam('', 4, bg[i]); {background slope}
      SetFunctionParam('', 5, os[i]); {background offset}
               
      PlotFunction(xFirst 1, xLast l[i] + 0.5, xStep 0, curveThickness 0.5 , curveRed 0x0000, curveGreen 0x0000, curveBlue 0x0000,
        legendEntry false, fillPattern 2);
   
      SetAxisProperties(axis 'X1', first 0, last pt);
  
      MoveTo(prleft + 98,prtop + 4 + (k-1)*30 + c + rectH - (CallFunction('', 1)-cbottom)/(ctop-cbottom)*rectH);
    	  SetTextStyle('Helvetica Regular',8,plain);
      NumberToString(os[i],s,false,-2);
    	  DrawText(s,0,false);
        
    end;
     
  EF := 1 + GetResult(fittedParameter, maxsamples*2 + 1);
  MoveTo(prright-140, 20 + nc * 30);
  NumberToString(EF,s,false,-4);
  SetTextStyle('Helvetica Regular',12,bold);
  DrawText('EA: ' + s , 0, false);
    
     
        
  end;