Lorenzian peak function

Dev Shed Forums Sponsor:

March 6th, 2012, 11:53 AM

Blackend

Registered User

Join Date: Mar 2012

Posts: 7

Time spent in forums: 1 h 21 m 45 sec
Reputation Power: 0

Lorenzian peak function

Hello!
I have got problem. I have to write a program that does the fit by Lorenz.

I'm stuck at the very beginning. I have a measurement file (. Txt). The file consists of two columns (x, y). I have to enter into the program (via button). Probably I have to use a an array or something similar. But I do not know how.
Then I want to get one after the other x, substitute into formulas. Then do the same with Y. And at the end insert the calculated x and y in the graph.
May I ask for help?

March 7th, 2012, 12:50 PM

Dimixx

Contributing User

Join Date: Jan 2010

Posts: 68

Time spent in forums: 23 h 53 m 47 sec
Reputation Power: 4

Hi
about loading if it's a *.txt file you can use a TStringList's method LoadFromFile.

Code:

var S : TStringList; 
begin s := TStringList.Create; 
s.LoadFromFile('YourFileName'); 
Memo1.Text := S.Text; 
end;

You can also use the Strings property of TStringList and then loop through the lines of the file. I'm don't know what's the formula you're trying do so if you have more questions I'll be glad to help you if I can.

Dimiter

March 7th, 2012, 01:07 PM

clivew

Contributing User

Join Date: Jan 2006

Location: Carlsbad, CA

Posts: 2,045

Time spent in forums: 1 Week 6 Days 2 h 37 m
Reputation Power: 382

To follow up on Dimixx suggestion, if you post a few lines from the text file that
illustrate how the x and y values on each line are formatted and separated I am sure
I or another member can take you a step or two further.

March 7th, 2012, 01:17 PM

Blackend

Registered User

Join Date: Mar 2012

Posts: 7

Time spent in forums: 1 h 21 m 45 sec
Reputation Power: 0

Thank for reply!

My code looks like this:

Code:

unit Pomiarowe;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, Styczne, ExtCtrls, TeeProcs, TeEngine, Chart, StdCtrls, Series, Symulacja,
  Spin;

type
  TForm4 = class(TForm)
    OpenDialog1: TOpenDialog;
    Chart1: TChart;
    Button1: TButton;
    Series1: TPointSeries;
    Button2: TButton;
    GroupBox1: TGroupBox;
    Panel1: TPanel;
    Edit1: TEdit;
    Panel2: TPanel;
    Edit2: TEdit;
    Panel3: TPanel;
    Edit3: TEdit;
    Panel7: TPanel;
    Edit7: TEdit;
    GroupBox2: TGroupBox;
    Panel5: TPanel;
    Edit5: TEdit;
    Panel6: TPanel;
    Edit6: TEdit;
    Button3: TButton;
    Label2: TLabel;
    Label1: TLabel;
    Series2: TPointSeries;
    Label4: TLabel;

    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
    u,v,w,uStart,vStart,wStart:Real;
    n,i:integer;
    X,Y: Real;
    function funk2(v:Real):Real;
    function funk3(w:Real):Real;
    function f1:Real;
    function f2(v:Real):Real;
    function f3(w:Real):Real;
    function Li(i:integer;u,v,w:Real):Real;
    function Lip(i:integer;v,w:Real):Real;
    procedure Fit;
  end;

var
  Form4: TForm4;

implementation
{$R *.DFM}


procedure TForm4.Button1Click(Sender: TObject);
var
  Dane: TextFile;
  x,y: Real;
begin
  with Form4 do
    begin
      if OpenDialog1.Execute then
        begin
          Series1.Clear;
          AssignFile(Dane, OpenDialog1.FileName);

          try
            Reset(Dane);

            while not EoF(Dane) do
              begin
                Readln(Dane, x, y);
                Form4.Series1.AddXY(x, y);
              end;
          finally
          CloseFile(Dane);
          end;
        end;
    end;
end;




function TForm4.Li(i:integer;u,v,w:Real):Real;
begin
 Result:=u*w/pi/(sqr(form4.x-v)+sqr(w));
end;

function TForm4.Lip(i:integer;v,w:Real):Real;
begin
 Result:=w/pi/(sqr(form4.x-v)+sqr(w));
end;

function TForm4.f1:Real;
var i:integer; licz,mian:Real;
begin
 Result:=0;
 licz:=0; mian:=0;
 for i:=0 to n do
 begin
  licz:= licz+lip(i,v,w)*form4.Y/abs(form4.Y);    
  mian:=mian+sqr(Lip(i,v,w))/abs(form4.Y);
     Result:=licz/mian;
 end;
end;

function TForm4.f2(v:Real):Real;
var i:integer;
begin
 Result:=0;
 for i:=0 to n do
 begin
 Result:=Result+(Li(i,u,v,w)-form4.Y)
 /abs(form4.Y)*sqr(Li(i,u,v,w))*(form4.X-v);
 end;
end;

function TForm4.f3(w:Real):Real;
var i:integer;
begin
 Result:=0;
 for i:=0 to n do
begin
 Result:=Result+(Li(i,u,v,w)-form4.Y)
 /abs(form4.Y)*sqr(Li(i,u,v,w));
 end;
end;

function TForm4.funk2(v:Real):Real;
begin
  Result:=0;
  if row_stycz(funk3,f3,wStart,w)then 
  Result:=f2(v);
  wStart:=w;
end;

function TForm4.funk3(w:Real):Real;
begin
 u:=f1;
 Result:=f3(w);
 uStart:=u;
end;




procedure TForm4.Fit;
var
        i:integer; chi2,r1,r2:Real;
begin
  if row_stycz(funk2,f2,vStart,v)then
  begin
   Edit1.text:=FloatToStr(u);
   Edit2.text:=FloatToStr(v);
   Edit3.text:=FloatToStr(w);
  Series1.Clear;
  Series2.Clear;
   chi2:=0;
   for i:=0 to n do
   begin

    r1:=form4.Y; r2:=(2*u/pi)*(w/(4*sqr(form4.X-v)+sqr(w)));
    Series1.AddXY(form4.X,r1);
    Series2.AddXY(form4.X,r2);


   chi2:=chi2+sqr(r2-r1)/abs(r1);

   end;
   chi2:=chi2/(n-3);
   Edit7.Text:=FloatToStr(chi2);
  end else ShowMessage('no solution');
end;




procedure TForm4.Button2Click(Sender: TObject);
begin
 vStart:=StrToFloat(Edit5.Text);
 wStart:=StrToFloat(Edit6.Text);
Fit;
end;



end.

I must change "procedure TForm4.Button1Click" for open txt file (via TStringList). But I want to do this via a button, and search the text file on hard disk. (I've got some problem with this, I am really beginner.)
After loading a text file, the program must take X and insert in place "form4.X" in every function. And the same to do with Y "form4.Y".
Finally, I want to display a text file on the chart (Series 1 - just data from *.txt; Series 2 - X, Y from fitting).

My text file looks like this:

Code:

28	3830.42023
30	3830.42023
33	3830.42023
35	3830.42023
38	3867.60878
....

Thanks for your help

March 7th, 2012, 03:31 PM

clivew

Contributing User

Join Date: Jan 2006

Location: Carlsbad, CA

Posts: 2,045

Time spent in forums: 1 Week 6 Days 2 h 37 m
Reputation Power: 382

OK.

A lot of complexity here and a number of Delphi "things" that could do with attention eventually.

However, starting with your specific issues

I have some questions.

1. Is your current Button1Click code actually filling Series1 with the correct x,y values?
If it is, then forget TStringLists etc. If not, then what is happening?

2. Why do you need/use TForm4.X and TForm4.Y if you already have the values
stored in series1? Of course I may have misread your intentions.
Would you not use Series1.XValue[n] and series1.YValue[n] in your loop?

Other quick comments.

You do not need and should not reference fields of TForm4 inside TForm4 with the Form4 or TForm4 prefix.
Most especially, do not reference them as Form4 rather than TForm4 as Form4 is
a specifically instantiated object of type TForm4 and may not be what
you want to address at all.
TForm4 is normally not needed but harmless (you could also use Self) Form4 could be dangerous.
One time you might include TForm4 would be to avoid ambiguity if,
for example, you wanted to be clear that you were referencing
TForm4.Caption and not Button1.Caption.
To take your button1Click code, the following would suffice:

Code:

procedure TForm4.Button1Click(Sender: TObject);
var
  Dane: TextFile;
  x,y: Real;
begin
      // Removed with Form4
      if OpenDialog1.Execute then
        begin
          Series1.Clear;
          AssignFile(Dane, OpenDialog1.FileName);
          try
            Reset(Dane);

            while not EoF(Dane) do
              begin
                Readln(Dane, x, y);
                Series1.AddXY(x, y);  // - Form4 removed here
              end;
          finally
            CloseFile(Dane);
          end;
        end;
end;

Also if any of the following are only used internally to TForm4 then make their declarations private.
u,v,w,uStart,vStart,wStart:Real;
n,i:integer;
X,Y: Real;

Finally (for now) it would probably be better to pass n into your procedures rather
than relying on it being set correctly at the global level.

Hope something in the above is helpful!

Clive

March 8th, 2012, 11:10 AM

Blackend

Registered User

Join Date: Mar 2012

Posts: 7

Time spent in forums: 1 h 21 m 45 sec
Reputation Power: 0

One more time thanks for reply

At the beginning I will answer questions.

1.
My current Button1Click code works fine. Chart fills with data (x,y). With your improvements also works correctly.

2.
I have to use a TForm4.X, TForm4.Y because the data from the file I need to perform calculations.
Probably I must use a loop that collects individual x and y. But unfortunately I do not know how to do it. I tried some tricks but it does not work. programming gives me a headache, but I have to do it.

Yes, I realize that should not be used TForm4 inside TForm4, etc.

I changed the declaration to private.

I changed everywhere form4.Y/form4.X to Y[n]/X[n] but the on the Series2 is only the last point. Though the calculations took only the last line of text file.

March 8th, 2012, 11:35 PM

clivew

Contributing User

Join Date: Jan 2006

Location: Carlsbad, CA

Posts: 2,045

Time spent in forums: 1 Week 6 Days 2 h 37 m
Reputation Power: 382

It is very hard to follow your code as it is full of math formulas that do not have their full context.

Leaving that aside for now:
1. I suggested accessing the values you wanted by using the XValue property as in
Series1.XValue[i] where i is the row of values you want.

2. I do not see where you even use the values you load.
As far as I can tell, you load the values in Button1Click but never use them.

Button2Click loads a couple of values and calls Fit.
Fit calls a procedure for which I see no code and, after that calls Series1.Clear which
immediately wipes out all the values you just loaded in Button1Click.

Another suggestion, if my math interpretation is correct.
I am using your original incorrect code for clarity of example.

If n can be a large number in function f1

licz:= licz+lip(i,v,w)*form4.Y/abs(form4.Y);

Would be faster (and clearer) if written as

Code:

if form4.Y > 0 then
  licz:= licz+lip(i,v,w)
else
  licz:= licz-lip(i,v,w);

I might be mistaken but isn't
form4.Y/abs(form4.Y)
simply going to evaluate to either 1 or -1
(unless form4.Y = 0 in which case you will get a divide by zero exception
which you are not catching).

Clive

March 9th, 2012, 01:34 AM

Blackend

Registered User

Join Date: Mar 2012

Posts: 7

Time spent in forums: 1 h 21 m 45 sec
Reputation Power: 0

I changed the code:

Code:

unit Pomiarowe;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, Styczne, ExtCtrls, TeeProcs, TeEngine, Chart, StdCtrls, Series, Symulacja,
  Spin;
type
  Wektor=array[0..1000] of real;
type
  TForm4 = class(TForm)
    OpenDialog1: TOpenDialog;
    Chart1: TChart;
    Button1: TButton;
    Series1: TPointSeries;
    Button2: TButton;
    GroupBox1: TGroupBox;
    Panel1: TPanel;
    Edit1: TEdit;
    Panel2: TPanel;
    Edit2: TEdit;
    Panel3: TPanel;
    Edit3: TEdit;
    Panel7: TPanel;
    Edit7: TEdit;
    GroupBox2: TGroupBox;
    Panel5: TPanel;
    Edit5: TEdit;
    Panel6: TPanel;
    Edit6: TEdit;
    Button3: TButton;
    Label2: TLabel;
    Label1: TLabel;
    Label4: TLabel;
    Series2: TPointSeries;

    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);

  private
    { Private declarations }
    u,v,w,uStart,vStart,wStart:Real;
    n,i:integer;
    x,y: wektor;
    //x,y: real;
  public
    { Public declarations }

    function funk2(v:Real):Real;
    function funk3(w:Real):Real;
    function f1:Real;
    function f2(v:Real):Real;
    function f3(w:Real):Real;
    function Li(i:integer;u,v,w:Real):Real;
    function Lip(i:integer;v,w:Real):Real;
    procedure Fit;
  end;

var
  Form4: TForm4;

implementation
{$R *.DFM}



procedure TForm4.Button1Click(Sender: TObject);
var
  Dane: TextFile;

begin
    begin
      if OpenDialog1.Execute then
        begin
          Series1.Clear;
          AssignFile(Dane, OpenDialog1.FileName);

          try
            Reset(Dane);

            while not EoF(Dane) do
              begin
                Readln(Dane, x[i], y[i]);
                Series1.AddXY(x[i], y[i]);
              end;
          finally
          CloseFile(Dane);
          end;
        end;
        end;
end;




function TForm4.Li(i:integer;u,v,w:Real):Real;
begin
 Result:=u*w/pi/(sqr(Series1.XValue[i]-v)+sqr(w));
end;

function TForm4.Lip(i:integer;v,w:Real):Real;
begin
 Result:=w/pi/(sqr(Series1.XValue[i]-v)+sqr(w));
end;

function TForm4.f1:Real;
var i:integer; licz,mian:Real;
begin
 Result:=0;
 licz:=0; mian:=0;
 for i:=0 to n do
 begin

 if series1.YValue[i] > 0 then
  licz:= licz+lip(i,v,w)
else
  licz:= licz-lip(i,v,w);

  mian:=mian+sqr(Lip(i,v,w))/abs(series1.YValue[i]);
     Result:=licz/mian;
 end;
end;

function TForm4.f2(v:Real):Real;
var i:integer;
begin
 Result:=0;
 for i:=0 to n do
 begin
 Result:=Result+(Li(i,u,v,w)-series1.YValue[i])/abs(series1.YValue[i])*sqr(Li(i,u,v,w))*(Series1.XValue[i]-v);
 end;
end;

function TForm4.f3(w:Real):Real;
var i:integer;
begin
 Result:=0;
 for i:=0 to n do
begin
 Result:=Result+(Li(i,u,v,w)-series1.YValue[i])/abs(series1.YValue[i])*sqr(Li(i,u,v,w));
 end;
end;

function TForm4.funk2(v:Real):Real;
begin
  Result:=0;
  if row_stycz(funk3,f3,wStart,w)then 
  Result:=f2(v);
  wStart:=w;
end;

function TForm4.funk3(w:Real):Real;
begin
 u:=f1;
 Result:=f3(w);
 uStart:=u;
end;




procedure TForm4.Fit;
var
        i:integer; chi2,r1,r2:Real;
begin
  if row_stycz(funk2,f2,vStart,v)then

  begin
   Edit1.text:=FloatToStr(u);
   Edit2.text:=FloatToStr(v);
   Edit3.text:=FloatToStr(w);

   chi2:=0;
   for i:=0 to n do
   begin

    r1:=Series1.YValue[i]; r2:=2*u/pi*(w/(4*sqr(Series1.XValue[i]-v)+sqr(w)));


    Series2.AddXY(Series1.XValue[i],r2);


   chi2:=chi2+sqr(r2-r1)/abs(r1);

   end;
   chi2:=chi2/(n-3);
   Edit7.Text:=FloatToStr(chi2);
  end else ShowMessage('no solution');
end;



procedure TForm4.Button2Click(Sender: TObject);
begin
 vStart:=StrToFloat(Edit5.Text);
 wStart:=StrToFloat(Edit6.Text);
Fit;
end;


end.

I also added an array, because probably I had to, right?
But still something is wrong. I still get only one point in the series2 - the first point, and that's all.
I put screenshot:
http://img402.imageshack.us/img402/7784/schowek01pr.jpg

March 9th, 2012, 11:28 AM

Dimixx

Contributing User

Join Date: Jan 2010

Posts: 68

Time spent in forums: 23 h 53 m 47 sec
Reputation Power: 4

Hi
I'm not very familiar with this type of Objects but it seem's to me that when you're adding your x and y in your Button1Click procedure
here :
while not EoF(Dane) do
begin
Readln(Dane, x[i], y[i]);
Series1.AddXY(x[i], y[i]);
end;
you don't change your index i.

Dimiter

#10

March 9th, 2012, 12:00 PM

clivew

Contributing User

Join Date: Jan 2006

Location: Carlsbad, CA

Posts: 2,045

Time spent in forums: 1 Week 6 Days 2 h 37 m
Reputation Power: 382

Deleted

Last edited by clivew : March 9th, 2012 at 12:03 PM. Reason: Responded before I saw the other new posts. My reply was out of date.

#11

March 9th, 2012, 12:22 PM

clivew

Contributing User

Join Date: Jan 2006

Location: Carlsbad, CA

Posts: 2,045

Time spent in forums: 1 Week 6 Days 2 h 37 m
Reputation Power: 382

Besides the error pointed out by Dimiter:

1. Where do you set the value of n? Did I miss it?
2. Do you ever initialize your two arrays with default values?
3. Do you check n is <= (rowsInFile-1) ?
4. Do you check that n is <= 1001?
5. Do you check that n is <= valid values in arrays?
(If something goes wrong with loading then (rowsInFile-1) may not
equal valid values in arrays).

Just comments to make your code more error resistant.

Clive

#12

March 9th, 2012, 01:11 PM

Blackend

Registered User

Join Date: Mar 2012

Posts: 7

Time spent in forums: 1 h 21 m 45 sec
Reputation Power: 0

1. Where do you set the value of n? Did I miss it?
I left the n value.

Now, I added a field to set the n value.

2. Do you ever initialize your two arrays with default values?
Yes.

3. Do you check n is <= (rowsInFile-1) ?
Yes. My n = 213.

4. Do you check that n is <= 1001?
Yes

5. Do you check that n is <= valid values in arrays?
Yes.

After setting the n values I received straight line in the Series2. So still there is something wrong, but there is little progress.

#13

March 10th, 2012, 09:22 AM

Blackend

Registered User

Join Date: Mar 2012

Posts: 7

Time spent in forums: 1 h 21 m 45 sec
Reputation Power: 0

Ok. Everything works well now!
I would like to thank you for all the tips and hints.

But I have one small problem. How to automate n values?
Now I made a window to manually set the value. But I would like to do it automatically.

Do you have any hints or ideas?

#14

March 10th, 2012, 10:26 AM

clivew

Contributing User

Join Date: Jan 2006

Location: Carlsbad, CA

Posts: 2,045

Time spent in forums: 1 Week 6 Days 2 h 37 m
Reputation Power: 382

Assuming n represents the number of rows in the file (a.k.a. the number of values in the series, one option would be
after loading the file into series1 to do

Code:

n := Series1.XValue.Count-1;

WARNING: I am not by a computer with Delphi installed. I have not confirmed that
the syntax and naming above is correct. However the solution to your question is.

Clive

#15

March 10th, 2012, 10:41 AM

Blackend

Registered User

Join Date: Mar 2012

Posts: 7

Time spent in forums: 1 h 21 m 45 sec
Reputation Power: 0

n := Series1.Count-1;
And everything works perfect

Again, many thanks to you clivew and Dimixx.

Lorenzian peak function

Developer Shed Advertisers and Affiliates