#include double f1(int id,double t,double *x[],double *v[]); double f2(int id,double t,double *x[],double *v[],double *B[]); //int imax=5E5;double endKickPar=0.5;//check kick int imax=6E5;double endKickPar=0.5;//check store double TK; int NFP=500*2; int NFP2=1000; int NF,IDUM,NFK,NFM; double FLR[500*2],FLZ[500*2],FLCRNT[500*2]; double FLRK[700],FLZK[700],FLCRNTK[700]; double FLRM[700],FLZM[700],FLCRNTM[700]; double Ikicker[50],posXc[50],posYc[50],posZc[50];//KICKER int cnKick;//kicker coil numbers double pi=TMath::Pi(); //--------------------------------------------------------------------- int cep12d_n10(double rk,double i,double &ak,double &ae,double &ill){ /* data a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10/ o 1.38629436111989d0, 0.0965735902811690d0, a 0.0308851465246711d0, 0.0149380448916805d0, b 8.79078273952743d-3, 6.18901033637687d-3, c 6.87489687449949d-3, 9.85821379021226d-3, d 7.97404013220415d-3, 2.28025724005875d-3, e 1.37982864606273d-4 / c data b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10/ * 0.5000d0 , 0.124999999999870d0, a 0.0703124996963957d0, 0.0488280347570998d0, b 0.0373774314173823d0, 0.0301204715227604d0, c 0.0239089602715924d0, 0.0154850516649762d0, d 5.94058303753167d-3, 9.14184723865917d-4, e 2.94078955048598d-5/ data ea0,ea1,ea2,ea3,ea4,ea5,ea6,ea7,ea8,ea9,ea10/ a 1.0000d0 , 0.443147180560990d0, b 0.0568051945617860d0, 0.0218317996015557d0, c 0.0115688436810574d0, 7.58395289413514d-3, d 7.77395492516787d-3, 0.0107350949056076d0, e 8.68786816565889d-3, 2.50888492163602d-3, f 1.53552577301013d-4/ c data eb0,eb1,eb2,eb3,eb4,eb5,eb6,eb7,eb8,eb9,eb10/ a 0.0d0 , 0.249999999999888d0, b 0.0937499997197644d0, 0.0585936634471101d0, c 0.0427180926518931d0, 0.0334833904888224d0, d 0.0261769742454493d0, 0.0168862163993311d0, e 6.50609489976927d-3, 1.00962792679356d-3, f 3.27954898576485d-5/ */ int ans=0; //data a0,a1,a2,a3,a4,a5,a6,a7,a8,a9,a10/ double a0=1.38629436111989e0; double a1=0.0965735902811690e0; double a2=0.0308851465246711e0; double a3=0.0149380448916805e0; double a4=8.79078273952743e-3; double a5=6.18901033637687e-3; double a6=6.87489687449949e-3; double a7=9.85821379021226e-3; double a8=7.97404013220415e-3; double a9=2.28025724005875e-3; double a10=1.37982864606273e-4; // data b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10/ double b0=0.5000e0; double b1=0.124999999999870e0; double b2=0.0703124996963957e0; double b3=0.0488280347570998e0; double b4=0.0373774314173823e0; double b5=0.0301204715227604e0; double b6=0.0239089602715924e0; double b7=0.0154850516649762e0; double b8=5.94058303753167e-3; double b9=9.14184723865917e-4; double b10=2.94078955048598e-5; // data ea0,ea1,ea2,ea3,ea4,ea5,ea6,ea7,ea8,ea9,ea10/ double ea0=1.0000e0; double ea1=0.443147180560990e0; double ea2=0.0568051945617860e0; double ea3=0.0218317996015557e0; double ea4=0.0115688436810574e0; double ea5=7.58395289413514e-3; double ea6=7.77395492516787e-3; double ea7=0.0107350949056076e0; double ea8=8.68786816565889e-3; double ea9=2.50888492163602e-3; double ea10=1.53552577301013e-4; // data eb0,eb1,eb2,eb3,eb4,eb5,eb6,eb7,eb8,eb9,eb10/ double eb0=0.0e0; double eb1=0.249999999999888e0; double eb2=0.0937499997197644e0; double eb3=0.0585936634471101e0; double eb4=0.0427180926518931e0; double eb5=0.0334833904888224e0; double eb6=0.0261769742454493e0; double eb7=0.0168862163993311e0; double eb8=6.50609489976927e-3; double eb9=1.00962792679356e-3; double eb10=3.27954898576485e-5; i=1; ill=0; double xm1,xm2,xm3,xm4,xm5,xm6,xm7,xm8,xm9,xm10,dalxm,alxm; double bzz; //if(rk>=0.0e0 && rk<=1.0e0) go to 5 if(rk<0.0e0 || rk>1.0e0){ ill=1; ak=0.0; ae=0.0; }else{ xm1=1.0e0-rk; xm2=xm1*xm1; xm3=xm2*xm1; xm4=xm3*xm1; xm5=xm4*xm1; xm6=xm5*xm1; xm7=xm6*xm1; xm8=xm7*xm1; xm9=xm8*xm1; xm10=xm9*xm1; dalxm=1.0e0/(double)xm1; alxm=log(dalxm); // bzz=b0 +b1*xm1+b2*xm2+b3*xm3+b4*xm4+b5*xm5+b6*xm6+b7*xm7+b8*xm8+b9*xm9+b10*xm10; ak= a0 +a1*xm1+a2*xm2+a3*xm3+a4*xm4+a5*xm5+a6*xm6+a7*xm7+a8*xm8+a9*xm9+a10*xm10 +bzz*alxm; bzz=eb0+eb1*xm1+eb2*xm2+eb3*xm3+eb4*xm4+eb5*xm5+eb6*xm6+eb7*xm7+eb8*xm8+eb9*xm9+eb10*xm10; ae=ea0 +ea1*xm1+ea2*xm2+ea3*xm3+ea4*xm4+ea5*xm5+ea6*xm6+ea7*xm7+ea8*xm8+ea9*xm9+ea10*xm10+bzz*alxm; ans=1; } return ans; } //--------------------------------------------------------------------- int cep12d(double RK,double I,double &AK,double &AE,double &ILL){ /* DATA AZ,A1,A2,A3,A4/ 1.38629436112D0,0.09666344259D0, A 0.03590092383D0,0.03742563713D0, B 0.01451196212D0 / DATA BZ,B1,B2,B3,B4/ 0.5D0 ,0.12498593597D0, A 0.06880248576D0,0.03328355346D0, B 0.00441787012D0/ DATA EAZ,EA1,EA2,EA3,EA4 / A 1.0D0 ,0.44325141463D0, B 0.06260601220D0,0.04757383546D0, C 0.01736506451D0/ DATA EBZ,EB1,EB2,EB3,EB4 / A 0.0D0 ,0.24998368310D0, B 0.09200180037D0,0.04069697526D0, C 0.00526449639D0/ */ int ans=0; double AZ=1.38629436112E0; double A1=0.09666344259E0; double A2=0.03590092383E0; double A3=0.03742563713E0; double A4=0.01451196212E0; double BZ=0.5E0; double B1=0.12498593597E0; double B2=0.06880248576E0; double B3=0.03328355346E0; double B4=0.00441787012E0; double EAZ=1.0E0; double EA1=0.44325141463E0; double EA2=0.06260601220E0; double EA3=0.04757383546E0; double EA4=0.01736506451E0; double EBZ=0.0E0; double EB1=0.24998368310E0; double EB2=0.09200180037E0; double EB3=0.04069697526E0; double EB4=0.00526449639E0; I=1; ILL=0; double XM1,XM2,XM3,XM4,DALXM,ALXM; double BZZ; //IF(RK>=0.0E0 && RK<=1.0E0) GO TO 5 if(RK<0.0E0 || RK>1.0E0){ ILL=1; AK=0.0; AE=0.0; }else{ XM1=1.0E0-RK; XM2=XM1*XM1; XM3=XM2*XM1; XM4=XM3*XM1; DALXM=1.0E0/(double)XM1; ALXM=log(DALXM); // BZZ=BZ+B1*XM1+B2*XM2+B3*XM3+B4*XM4; AK= AZ+A1*XM1+A2*XM2+A3*XM3+A4*XM4 + BZZ*ALXM; // BZZ=EBZ+EB1*XM1+EB2*XM2+EB3*XM3+EB4*XM4; AE=EAZ+EA1*XM1+EA2*XM2+EA3*XM3+EA4*XM4+BZZ*ALXM; ans=1; } return ans; } //--------------------------------------------------------- int bfield(double CR,double CZ,double CI,double RI,double ZJ,double &BR,double &BZ,double &APHI){ // COMMON/COND/ ILL // DOUBLE PRECISION RK, ELPK,ELPE // // (CR,CZ)--- POSITION OF COIL // (RI,ZJ)--- MEASUREMENT POSITION // (BR,BZ)--- MAGNETIC FIELD // (APHI) --- VECTOR POTENTIAL // RETURN int cep=0; double XMU=4.E-7; double S =RI*RI+CR*CR+(ZJ-CZ)*(ZJ-CZ); double P =RI*CR+RI*CR; double RK=(P+P)/(double)(S+P); double PSI; // printf("RK=%lf,S=%lf,P=%lf\n",RK,S,P); double ELPK,ELPE,ILL; //cep=cep12d(RK,1,ELPK,ELPE,ILL); cep=cep12d_n10(RK,1,ELPK,ELPE,ILL);//2022Apr25 // printf("cep-out: RK=%e ELPK=%e ELPE=%e ILL=%e\n",RK,ELPK,ELPE,ILL);getchar(); if(ILL==0){ RK=sqrt(RK); BZ =XMU*CI/(double)(2.E0* sqrt(S+P))*(ELPK-(S-2.E0*CR*CR)/(double)(S-P)*ELPE); BR =XMU*CI/(double)(2.E0* sqrt(S+P))*(ZJ-CZ)/(double)RI*(-ELPK+ S/(double)(S-P)*ELPE); PSI=CI/(double)RK* sqrt(RI*CR)*((1.E0-RK*RK/(double)2.E0)*ELPK-ELPE); APHI=XMU*PSI/(double)RI; // printf("ILL=%e,RK=%e,BZ=%e,BR=%e\n",ILL,RK,BZ,BR);getchar(); } //printf("ILL=%lf,RK=%lf,S=%lf,P=%lf\n",ILL,RK,S,P); //1000 FORMAT(1H ,15(1H*),' WARNING IN FILD ILL=',I5, 5(1H*),3(1PE10.3)) //<--? return cep; } //--------------------------------------------------------------- int bflfit(double RM,double ZM,double &BR,double &BZ,double &APHI){ int ans=0; int ans1=0; BR=0.0E0; BZ=0.0E0; APHI=0.0E0; double DDD,DBR,DBZ,DAPHI; for(int i=0;i420)getchar(); } // printf("IN----BR=%lf BZ=%lf\n",BR,BZ);getchar(); ans=1; return ans; } //--------------------------------------------------------------- int bflfit2(double RM,double ZM,double &BR,double &BZ,double &APHI,double I[],double posXc[],double posYc[],int cn){ int ans=0; BR=0.0E0; BZ=0.0E0; APHI=0.0E0; double DDD,DBR,DBZ,DAPHI; for(int i=0;i(NFM-9)){ FLCRNTM[i]=FLCRNTM[i]*0.48;// //FLCRNTM[i]=FLCRNTM[i]*0.;// } } fclose(infile); ans=1; return ans; } double f1(int id,double t,double x[3],double v[3]) { return v[id]; } double f2(int id,double t,double x[3],double v[3],double B[3]) { double ans[3]; ans[0]=v[1]*B[2]-v[2]*B[1]; ans[1]=v[2]*B[0]-v[0]*B[2]; ans[2]=v[0]*B[1]-v[1]*B[0]; return ans[id]; } //////////////KICK-SPATIAL int kickInput(double facB,int kickTyp){ FILE *infile; char sdmy[256]; int idmy; char iname[256]; if(kickTyp==9 )sprintf(iname,"KickCoil_file/Reverse_20220215_typ%d.txt",kickTyp); if(kickTyp==7 )sprintf(iname,"KickCoil_file/Reverse_20220202_typ%d.txt",kickTyp); if(kickTyp==11 || kickTyp==14)sprintf(iname,"Reverse_20220329_typ%d.txt",kickTyp); if(kickTyp>=15)sprintf(iname,"KickCoil_file/Reverse_20220412_typ%d.txt",kickTyp); if(kickTyp==23 || kickTyp==24)sprintf(iname,"KickCoil_file/Reverse_20220427_typ%d.txt",kickTyp); if(kickTyp>=25 && kickTyp<28)sprintf(iname,"KickCoil_file/Reverse_20220524_typ%d.txt",kickTyp); if(kickTyp>=28 && kickTyp<31)sprintf(iname,"KickCoil_file/Reverse_20220619_typ%d.txt",kickTyp); //sprintf(iname,"Reverse_20220412_typ%d.txt",kickTyp); infile=fopen(iname,"r"); //infile=fopen("Reverse_20220329_typ11.txt","r"); //fscanf(infile,"%s",sdmy);//printf("aaa\n");getchar(); int cc=0; for(int k=0;k<2;++k){ fgets(sdmy,256,infile);//printf("%s",sdmy);//getchar(); while(fscanf(infile,"%d %lf %lf %lf",&idmy,&posXc[cc],&posZc[cc],&Ikicker[cc])==4){ Ikicker[cc]=Ikicker[cc]*facB; cc++; } } fclose(infile); cnKick=cc; // printf("cnKick=%d\n",cnKick); return 0; } int kickSpatial3(double facB,double ratio,double x[3],double &BRkick0,double &BZkick0,int kickTyp){ int ff=0; int home=0; double ldmy; int idmy; // facB=1; kickInput(facB,kickTyp); // printf("cnKick=%d\n",cnKick); /////////////check B-field///////// double ZM,RM; double BR,BZ,APHI; ZM=x[2]; RM=sqrt(pow(x[0],2)+pow(x[1],2)); bflfit2(RM,ZM,BR,BZ,APHI,Ikicker,posXc,posZc,cnKick); BRkick0=BR*1E4;//Gauss BZkick0=BZ*1E4; return 1; } //////////////KICK-TOOL int kickMag(double t,double fac,double x[3],double BRkick0,double BZkick0,double &BRkick,double &BZkick){ BRkick0=BRkick0*1E-4; BZkick0=BZkick0*1E-4; double theta=t/(double)(TK*1E-9)*2*pi; BRkick=BRkick0*sin(theta); BZkick=BZkick0*sin(theta); return 1; } ////Solenoid z-axis void MagFieldKick() { int kickTyp=11; int weakTyp=2; int now_i=0; int now_j=0; int ff=0; FILE *ofile,*infile; // FILE *ofile2; char oname2[256]; char oname[256]; char iname[256]; char sdmy[256]; double ldmy; int idmy; double xtmp,ytmp,ztmp,ts; double timeStop; double t=0; double x[3]; double B[3]={0,0,3}; double xp[128],yp[128],zp[128],tp[128]; double bx[128],by[128],bz[128],br[128]; /// INJECTION INPUT FILE double t0,dtfac; double facB,start_t0;//weak_0.45 mod-3 redo 20190018 New Kick2_3 sprintf(iname,"kickInp.txt"); infile=fopen(iname,"r"); fscanf(infile,"%s",sdmy);printf("1--->%s\n",sdmy); fscanf(infile,"%lf,%lf",&facB,&TK); fscanf(infile,"%s",sdmy);printf("2-->%s\n",sdmy); int evtCnt=0; while(fscanf(infile,"%lf,%lf,%lf,%lf",&xtmp,&ytmp,&ztmp,&ts)==4){ xp[evtCnt]=xtmp; //m yp[evtCnt]=ytmp; //m zp[evtCnt]=ztmp; //m tp[evtCnt]=ts*1E-9; //sec evtCnt++; if(evtCnt>127)break; } fclose(infile); printf("evtCnt=%d\n",evtCnt); ////////////////////point_INPUT END///////////// ///////////////read FLDATA////WEAK FOCUS/////// double ratio=facB; double RM,ZM,BR,BZ,APHI; int read=0; int fitp=0; read=readfldata_20210727();//STATIC weak+main if(read!=1){ printf("FLDATA read err\n"); exit(0); } /////////////set kickB-field///////// double end_kick=0.5*TK*1E-9;//sec //halfsine double BRkick,BZkick; double BRkick0,BZkick0; double kickTime; double int_time; int kick_status; for(int i=0;i=end_kick)kick_status=2; /////////////Get B-WEAK-field///////// fitp=bflfitM(RM,ZM,BR,BZ,APHI);//CALC-ABE-SAN WEAK FIELD 20210727 //printf("BY=%lf BR=%lf\n",BY,BR);getchar(); /////////////SUM B-KICK+B-WEAK-field///////// br[i]=(BR+BRkick); //br[i]=(BRkick); //kicker only bx[i]=(BR+BRkick)*x[0]/(double)RM; by[i]=(BR+BRkick)*x[1]/(double)RM; bz[i]=BZ+BZkick; }//END for(evtCnt) ///////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); TGraph *gr_BRall; gr_BRall = new TGraph(evtCnt,br,zp); TCanvas *c = new TCanvas("c","br-z",200,100,500,500); c->SetLeftMargin(0.13); //TH2F* frame = new TH2F("frame"," ",10, -0.4, 0.4, 100, 0., 0.9); TH2F* frame = new TH2F("frame"," ",10, -0.1, 0.05, 100, 0., 0.9); frame->SetTitle(""); frame->GetXaxis()->SetLabelSize(0.04); frame->GetXaxis()->SetTitleSize(0.05); frame->GetYaxis()->SetLabelSize(0.04); frame->GetYaxis()->SetTitleSize(0.05); frame->GetXaxis()->SetTitle("br[T]"); frame->GetYaxis()->SetTitle("z[m]"); frame->GetYaxis()->SetTitleOffset(1.2); frame->SetStats(0); frame->Draw(); gr_BRall->SetMarkerStyle(20); gr_BRall->SetMarkerColor(1); gr_BRall->SetLineWidth(2); gr_BRall->Draw("pl"); c->SetGrid(); }