Subversion Repositories f9daq

#include <stdlib.h>

#include <stdio.h>

#include <TString.h>

#include <TGraph.h>

#include <TCanvas.h>

#include <TMultiGraph.h>

#include <math.h>

int IUdraw(const char *datoteka =  "C:/home/data/meritev-iv.dat"){

auto c1 = new TCanvas("c1","A Simple Graph Example",200,10,1400,1000);

double sigmaGauss = 0.3;

double xError = 0.5;

double rangeX = 5;

//c1->SetLogy(); 

//auto mg = new TMultiGraph();

// printf("#%s#\n", datoteka);

 /*

string s = datoteka;

s.erase(s.find_last_of("."), string::npos);

fprintf(stderr,"%s\n", s.c_str());

std::string key ("/");

string ime;

std::size_t found = s.rfind(key);

if (found!=std::string::npos) {

  ime = s.substr(found+1,s.length());

}

*/

TString ime = gSystem->BaseName(datoteka);

TString path = gSystem->DirName(datoteka);

ime.Remove(ime.Index(".dat"));

fprintf(stderr,"Ime=%s\n", ime.Data());

TString formatted = path + TString("/") + ime; // + TString(".png");

//formatted.Form("%s.png", );

fprintf(stderr,"png=%s\n", formatted.Data());

TString formattedGraph1;

TString formattedGraph2;

TString formattedGraph3;

TString formattedGraph4;

TString formattedGraph5;

TString formatted1;

TString formatted2;

TString formatted3;

TString formatted4;

TString formatted5;

for (Int_t d = 0; d<1; d++) {

  FILE * fp=fopen(datoteka,"r");

  if (!fp) continue;

  float f[4];

  Int_t j=0;  

  Int_t ndim=400;

  char line[400];  

  double fpXX[250];

  double fpYY[250];

//  Int_t k = 0;

  while (fgets(line,ndim,fp)!=NULL) {

        if (line[0] == '#') continue;

    sscanf(line,"%f%f%f%f",&f[0],&f[1],&f[2],&f[3]);

//      printf("%f %f %f %f \n",f[0],f[1],f[2],f[3]);

        fpXX[j]=f[2];

        fpYY[j]=f[3];

        j++;

  }

  // float *fpx =    new float[j+1];

  // float *fpy =    new float[j+1];

  // for (Int_t i=0; i<j-3; i++) { 

    // if (fpXX[i]!= 0) {

                // k++;

                // fpx[i]=fpXX[i+2];

                // fpy[i]=fpYY[i+2];

                // printf("%2.6f\t",fpx[i]);

                // printf("\n");

                // printf("%2.6f\t",fpy[i]);

                // printf("\n");

        // }

  // }

  fclose(fp);

  auto graf11 = new TGraph();

  graf11->SetMarkerStyle(20);

  graf11->SetDrawOption("AP");

  graf11->SetLineWidth(3);

  graf11->SetFillStyle(0);

  auto graf12 = new TGraph();

  graf12->SetMarkerStyle(20);

  graf12->SetDrawOption("AP");

  graf12->SetLineWidth(3);

  graf12->SetFillStyle(0);

  auto graf2 = new TGraph();

  graf2->SetMarkerStyle(20);

  graf2->SetDrawOption("AP");

  graf2->SetLineWidth(3);

  graf2->SetFillStyle(0);

  auto graf3 = new TGraph();

  graf3->SetMarkerStyle(20);

  graf3->SetDrawOption("AP");

  graf3->SetLineWidth(3);

  graf3->SetFillStyle(0);

  auto graf4 = new TGraph();

  graf4->SetMarkerStyle(20);

  graf4->SetDrawOption("AP");

  graf4->SetLineWidth(3);

  graf4->SetFillStyle(0);

  auto graf5 = new TGraph();

  graf5->SetMarkerStyle(20);

  graf5->SetDrawOption("AP");

  graf5->SetLineWidth(3);

  graf5->SetFillStyle(0);

  for (int i=0; i<j; i++) {

        if (fpYY[i]<=0) fpYY[i] = 1e-12;  

        graf11->SetPoint(i,fpXX[i],log(fpYY[i]));

        graf12->SetPoint(i,fpXX[i],log(fpYY[i]));

        graf5->SetPoint(i,log(fpYY[i]),fpXX[i]);

  }

  double max2 = 0;

  double max2x = 0;

  int k = 0;

  for (int i=0; i<j-1; i++) {

        if (fpYY[i]<=0) fpYY[i] = 1e-12;

        double x2 = (fpXX[i+1]+fpXX[i])/2;

        double y2 = (log(fpYY[i+1])-log(fpYY[i]))/(fpXX[i+1]-fpXX[i]);

        graf2->SetPoint(i,x2, y2);

        if (y2 > 1.1) {

                double y3 = 1 / y2;

                graf3->SetPoint(k,x2, y3);

                k++;

        }

        if (y2 > max2) {

                max2 = y2;

                max2x = x2;

        }

  }

  double max4 = 0;

  double max4x = 0;

  double x4stari = 0;

  double y4stari = 0;

  double x4novi = 0;

  double y4novi = 0;

  for (int i=0; i<j-2; i++) {

        if (fpYY[i]<=0) fpYY[i] = 1e-12;

        double x4 = (fpXX[i+2]+fpXX[i])/4+fpXX[i+1]/2;

        double y4 = ((log(fpYY[i+2])-log(fpYY[i+1]))/(fpXX[i+2]-fpXX[i+1]) - (log(fpYY[i+1])-log(fpYY[i]))/(fpXX[i+1]-fpXX[i]))/((fpXX[i+2]-fpXX[i])/2);

        if (i==0) {

                x4novi = x4;

                y4novi = y4;

        }

        else {

                x4novi = (x4 + x4stari)/2;

                y4novi = (y4 + y4stari)/2;

        }

        graf4->SetPoint(i,x4novi, y4novi);

        x4stari = x4;

        y4stari = y4;

        if (y4 > max4) {

                max4 = y4;

                max4x = x4;

        }

  }

}

TF1 *fa1 = new TF1("fa1","pol1",61,67);

fa1->SetParameters(-28.4467,0.0875386);

graf11->Fit("fa1", "R");

graf11->Draw("APL");

c1->Update();

TF1 *fa2 = new TF1("fa2","pol1",69.7,70.1);

graf11->Fit("fa2", "R+");

//graf12->Draw("APL");

printf("p0=%e, p1=%e\n", fa1->GetParameter(0), fa1->GetParameter(1));

double par0 = fa1->GetParameter(0);

double par1 = fa1->GetParameter(1);

printf("p2=%e, p3=%e\n", fa2->GetParameter(0), fa2->GetParameter(1));

double par2 = fa2->GetParameter(0);

double par3 = fa2->GetParameter(1);

double x1 = (par0 - par2)/(par3 - par1);

formattedGraph1.Form("I ( V ) -> %s, V(1) = %f ; V [V]; ln(I) [A]", ime.Data(), x1);

graf11->SetTitle(formattedGraph1);

graf12->SetTitle(formattedGraph1);

//graf13->Draw("APL");

formatted1.Form("%s-1.png", formatted.Data());

c1->SaveAs(formatted1);

auto c2 = new TCanvas("c2","A Simple Graph Example",200,10,1400,1000);

graf2->Draw("APL");

graf2->GetXaxis()->SetRangeUser(x1 - rangeX, x1 + rangeX);

TF1 *fa2 = new TF1("fa2","gaus(0)+pol1(3)",max2x - sigmaGauss, max2x + sigmaGauss);

 fa2->SetParameters(60,x1,0.1339,0,0);

 fa2->SetParLimits(1, max2x - sigmaGauss, max2x + sigmaGauss);

// fa1->SetParLimits(1, 10**(-14), 10**(-10));

 graf2->Fit("fa2", "R");

// printf("p20=%e, p21=%e, p22=%e, max2x=%f\n", fa2->GetParameter(0), fa2->GetParameter(1), fa2->GetParameter(2), max2x);

// double par2 = fa2->GetParameter(2);

// double par3 = fa2->GetParameter(3);

double x2 = fa2->GetParameter(1);

formattedGraph2.Form("I ( V ) -> %s, V(2) = %f ; Voltage [V]; d(ln(I)/dV)", ime.Data(), x2);

graf2->SetTitle(formattedGraph2);

formatted2.Form("%s-2.png", formatted.Data());

gStyle->SetOptFit(1);

c2->SaveAs(formatted2);

auto c3 = new TCanvas("c3","A Simple Graph Example",200,10,1400,1000);

graf3->Draw("APL");

graf3->Fit("pol3");

graf3->GetXaxis()->SetRangeUser(x1 - rangeX, x1 + rangeX);

TF1 *fa3 = new TF1("fa3", "pol1", x2, x2 + 0.5);

fa3->SetParameters(-30.6306,0.4398);

//fa3->SetParLimits(0, 10**(-14), 10**(-10));

//fa3->SetParLimits(1, 10**(-14), 10**(-10));

graf3->Fit("fa3", "R");

double n3 = fa3->GetParameter(0);

double k3 = fa3->GetParameter(1);

double x3 = - n3 / k3;

formattedGraph3.Form("I ( V ) -> %s, V(3) = %f ; Voltage [V]; 1/(d(ln(I)/dV))", ime.Data(), x3);

graf3->SetTitle(formattedGraph3);

//printf("p30=%e, p31=%e\n", fa3->GetParameter(0), fa3->GetParameter(1));

formatted3.Form("%s-3.png", formatted.Data());

c3->SaveAs(formatted3);

auto c4 = new TCanvas("c4","A Simple Graph Example",200,10,1400,1000);

graf4->Draw("APL");

graf4->GetXaxis()->SetRangeUser(x1 - rangeX, x1 + rangeX);

TF1 *fa4 = new TF1("fa4","gaus(0)",max4x - sigmaGauss, max4x + sigmaGauss);

 fa4->SetParameters(60,x2,0.1339);

 fa4->SetParLimits(1, x2 - xError, x2 + xError);

 fa4->SetParLimits(0, 0, 1000);

 graf4->Fit("fa4", "R");

// printf("p40=%e, p41=%e, p42=%e, max4x=%f\n", fa4->GetParameter(0), fa4->GetParameter(1), fa4->GetParameter(2), max4x);

// double par2 = fa2->GetParameter(2);

// double par3 = fa2->GetParameter(3);

double x4 = fa4->GetParameter(1);

formatted4.Form("%s-4.png", formatted.Data());

// TString stringx4;

// stringx4.Form("V(4) = %f", x4);

   // TText *t = new TText(60,70,stringx4);

   // t->SetTextAlign(22);

   // t->SetTextColor(1);

   // t->SetTextFont(43);

   // t->SetTextSize(20);

   // t->Draw();

formattedGraph4.Form("I ( V ) -> %s, V(4) = %f ; Voltage [V]; d^2(ln(I)/dV^2)", ime.Data(), x4);

graf4->SetTitle(formattedGraph4);

gPad->Modified(); gPad->Update();

c4->SaveAs(formatted4);

auto c5 = new TCanvas("c5","A Simple Graph Example",200,10,1400,1000);

graf5->Draw("APL");

TF1 *fa5 = new TF1("fa5","-[0] - [1]*x - [2]*x**2",par0 + par1*x1 +4,par0 + par1*x1+11);

// fa5->SetParameters(28.4467,0.0875386,4);

graf5->Fit("fa5", "R+");

printf("p4=%e, p5=%e, p6=%e\n", fa5->GetParameter(0), fa5->GetParameter(1), fa5->GetParameter(2));

double par4 = fa5->GetParameter(0);

double par5 = fa5->GetParameter(1);

double par6 = fa5->GetParameter(2);

double xx5 = - (double)(par5/(2*par6));

double x5 = - par4 - par5 * xx5 - par6 * xx5**2;

// double f5(double p);

// double f5(double p) {

        // double par0, par1, par4, par5, par6;

    // double a=-(double)(par0/par1) + (double)(1/par1) * p + par4 + par5 * p + par6 * p**2;

    // return a;

// }

// cout.precision(4);        

// cout.setf(ios::fixed);

// double a5,b5,cc5,e5,faa5,fb5,fc5; 

// a5 = x1-5;

// b5 = x1+5;

// e5 = 0.001;

    // if (f5(a5)*f5(b5)>0) print("No root exists between a and b.\n");

    // else {

    // while (fabs(a5-b5)>=e5) {

        // cc5=(a5+b5)/2.0;        //bisect the interval and find the value of c

        // faa5=f5(a5);

        // fb5=f5(b5);

        // fc5=f5(cc5);  

        // if (fc5==0) {

            // printf("The root of the equation is %f \n",cc5); 

            // break;

        // }

        // if (fa5*fc5>0) a5=cc5;

        // else if (fa5*fc5<0) b5=cc5;

    // }

    // }            //The loop ends when the difference between a and b becomes less than the desired accuracy ie now the value stored in 'c' can be called the approximate root of the equation         

    // printf("The root of the equation is %f \n",cc5);  

formattedGraph5.Form("V ( I ) -> %s, V(5) = %f ; ln(I) [A] ; V [V]", ime.Data(), x5);

graf5->SetTitle(formattedGraph5);

formatted5.Form("%s-5.png", formatted.Data());

c5->SaveAs(formatted5);

printf("RESULT: V(1) = %f, V(2) = %f, V(3) = %f, V(4) = %f, V(5) = %f", x1, x2, x3, x4, x5);

return 0;

}