int npoint=256*1;
void Bio10(){
 double mu0=1.256637062*1E-6;
 double x[100][npoint],y[100][npoint],z[100][npoint];//max 100 coils
 double Hx[100],Hy[100],Hz[100];
 double xm[100],ym[100],zm[100];//reference
 int nm=10;
 int jn=npoint;
 double cir_R=0.35;//[m]
 double I=-800;//[A]
 double pi=TMath::Pi();

 TPolyMarker3D *check;
  check=new TPolyMarker3D(256);
 TPolyLine3D *checkCoil[100];
  for(int k=0;k<100;++k)checkCoil[k]=new TPolyLine3D(256*40);
//reference point
 for(int i=0;i<nm*1;++i){
    xm[i]=0;
    ym[i]=0.;
    zm[i]=0.5-1.0*i/(double)nm;
      check->SetPoint(i,xm[i],ym[i],zm[i]);
 }
//ring1-100
int nn=0;
printf("how many coils do you want? (upto 50 ) Please put number=");scanf("%d",&nn);
if(nn<0 || nn>50){
	printf("wrong number!\n");
        exit(0);
}
nn=nn*2;
double dz=0.01;
 for(int k=0;k<nn;++k){
 for(int j=0;j<jn;++j){
    double theta=2*pi*j/(double)jn;
    x[k][j]=cir_R*cos(theta);
    y[k][j]=cir_R*sin(theta);
    if(k%2==0)z[k][j]=dz+dz*k;
    if(k%2==1)z[k][j]=-1*dz-1*dz*k;
      checkCoil[k]->SetPoint(j,x[k][j],y[k][j],z[k][j]);
 }
 }
//ring1-100
double Bioz=0;
I=I/(double)nn;
 for(int i=0;i<nm;++i){

   Hx[i]=0; Hy[i]=0; Hz[i]=0;
   for(int k=0;k<nn;++k){
    for(int j=0;j<(jn-1);++j){
      double ds[3];
      double rr[3];
      double dHx,dHy,dHz;

      ds[0]=x[k][j+1]-x[k][j];
      ds[1]=y[k][j+1]-y[k][j];
      ds[2]=z[k][j+1]-z[k][j];
      rr[0]=xm[i]-x[k][j];
      rr[1]=ym[i]-y[k][j];
      rr[2]=zm[i]-z[k][j];

      double rr0=sqrt(rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2]);
      dHx=I*(ds[1]*rr[2]-ds[2]*rr[1])/(double)pow(rr0,3);
      dHy=I*(ds[2]*rr[0]-ds[0]*rr[2])/(double)pow(rr0,3);
      dHz=I*(ds[0]*rr[1]-ds[1]*rr[0])/(double)pow(rr0,3);
//      printf("dHx=%e dHy=%e dHz=%e\n",dHx,dHy,dHz);getchar();
      Hx[i]=Hx[i]+dHx;
      Hy[i]=Hy[i]+dHy;
      Hz[i]=Hz[i]+dHz;
    }
   }
//Field flux

    Hx[i]=mu0*Hx[i]/(double)(4*pi);//tesla
    Hy[i]=mu0*Hy[i]/(double)(4*pi);
    Hz[i]=mu0*Hz[i]/(double)(4*pi);
    printf("point at %d x=%lf y=%lf z=%lf Bx=%e By=%e Bz=%e\n",i,xm[i],ym[i],zm[i],Hx[i],Hy[i],Hz[i]);
    if(sqrt(pow(zm[i],2))<1E-3)Bioz=Hz[i];
   } 

//Omajinai to write graphs
 TGraph *gr_Bx,*gr_By,*gr_Bz;
  gr_Bx=new TGraph(nm,zm,Hx);
  gr_By=new TGraph(nm,zm,Hy);
  gr_Bz=new TGraph(nm,zm,Hz);
//H=I/2r   circle
//H=I/2pir  line 

    double B_test=mu0*I/(double)(2*cir_R);
    printf("B_check=%e Biot-savart low=%e\n",B_test, Bioz);
    printf("ration=%e\n",(Bioz-B_test)/B_test);
///////Graph
gROOT->SetStyle("Plain");
gROOT->ForceStyle();
gStyle->SetOptStat(0);
gStyle->SetCanvasBorderMode(0);
gStyle->SetPadColor(0);
gStyle->SetTitleFillColor(0);
gStyle->SetStatColor(0);
//
 gStyle->SetCanvasBorderMode(0);//remove yellow line
 gStyle->SetCanvasColor(0);
 gStyle->SetFrameLineColor(1);
 gStyle->SetPadColor(0);
 gStyle->SetTitleFillColor(0);
 gStyle->SetStatColor(0);
 TCanvas *c2 = new TCanvas("c2","coil",200,100,1000,700);
 TH3F *framecc2;
   framecc2 = new TH3F("framecc2","",10, -0.4, 0.4,10, -0.4, 0.4, 10,-0.7,0.7);
         framecc2->GetXaxis()->SetLabelSize(0.04);
         framecc2->GetXaxis()->SetTitleSize(0.04);
         framecc2->GetYaxis()->SetLabelSize(0.04);
         framecc2->GetYaxis()->SetTitleSize(0.04);
         framecc2->GetZaxis()->SetLabelSize(0.04);
         framecc2->GetZaxis()->SetTitleSize(0.04);
         framecc2->GetXaxis()->SetTitle("x [m]");
         framecc2->GetYaxis()->SetTitle("y [m]");
         framecc2->GetZaxis()->SetTitle("z [m]");
         framecc2->GetXaxis()->SetTitleOffset(1.2);
         framecc2->GetYaxis()->SetTitleOffset(1.2);
         framecc2->GetZaxis()->SetTitleOffset(1.6);
              framecc2->SetStats(0);
              framecc2->Draw();

    c2->SetLeftMargin(0.15);
    c2->SetRightMargin(0.05);
    c2->SetBottomMargin(0.10);

   check->SetMarkerStyle(24);
   check->SetMarkerColor(2);
   check->Draw("p");
   for(int k=0;k<nn;++k){
     checkCoil[k]->SetLineColor(6);
     checkCoil[k]->Draw("l");
   }

   ////////new plot of Bx, By and Bz 20200416
 TCanvas *c3 = new TCanvas("c3","BB",200,100,1000,700);
 TH2F *framecc3;
   framecc3 = new TH2F("framecc3","",10, -0.7, 0.7,10, -2E-3,1E-3);
         framecc3->GetXaxis()->SetLabelSize(0.04);
         framecc3->GetXaxis()->SetTitleSize(0.04);
         framecc3->GetYaxis()->SetLabelSize(0.04);
         framecc3->GetYaxis()->SetTitleSize(0.04);
         framecc3->GetXaxis()->SetTitle("z [m]");
         framecc3->GetYaxis()->SetTitle("B [tesla]");
         framecc3->GetXaxis()->SetTitleOffset(1.2);
         framecc3->GetYaxis()->SetTitleOffset(1.5);
              framecc3->SetStats(0);
              framecc3->Draw();

    c3->SetLeftMargin(0.15);
    c3->SetRightMargin(0.05);
    c3->SetBottomMargin(0.10);
    
gr_Bx->SetMarkerStyle(20);
gr_By->SetMarkerStyle(20);
gr_Bz->SetMarkerStyle(20);
gr_Bx->SetMarkerColor(2);
gr_By->SetMarkerColor(3);
gr_Bz->SetMarkerColor(kViolet-1);
gr_Bx->SetLineColor(2);
gr_By->SetLineColor(3);
gr_Bz->SetLineColor(kViolet-1);
gr_Bx->Draw("pl");
gr_By->Draw("pl");
gr_Bz->Draw("pl");
c3->SetGrid();
   //
}