void currents_cool(void) { // // plot the steady state solution to the heat transer problem with a source from [-sigma,sigma] in a uniform cooling bath // // #include #include gROOT->Reset(); //TTree *Bfield = (TTree *) gROOT->FindObject("Bfield"); gStyle->SetPalette(1,0); gStyle->SetOptStat(kFALSE); gStyle->SetOptFit(kFALSE); // gStyle->SetOptFit(1111); gStyle->SetPadRightMargin(0.15); gStyle->SetPadLeftMargin(0.15); gStyle->SetPadBottomMargin(0.15); gStyle->SetFillColor(0); // char string[256]; char filename[80]; Int_t j,jj; #define npts 4; #define jmax 10001; #define kmax 4; # define pi 3.14159; // input data set first point (actually -100) to 100 so that it is not plotted. // use offset on x-axis to make data visible and give negative y-values for reference (not plotted) Double_t dummyx[npts]={0.1,2.,10.,5000}; Double_t dummyy[npts]={-10000,-10000,-10000,-10000}; Double_t xmin=0.1; Double_t xmax=2; //Double_t xmin=0.5; //Double_t xmax=10; Double_t ymin=-3000; Double_t ymax=10000; Double_t gamma0=0.34; // Double_t gamma0=0.5; Double_t gamma=gamma0; Double_t turns=4608; Double_t climit=10000; Double_t L1= 26. ; // inductance of solenoid, units are H minus one coil Double_t L2 =7.3e-6; // self-inductance of loop, units are H, computed from single loop equation Double_t L2 =L1/(turns*turns); // scale to inductance of full solenoid divided by the number of turns. Double_t L2=10e-6; // calculation by Eugene for the self-inductance. // Double_t L3 = L1/4600; // self-inductance of one coil due to solenoid Double_t L3 = 0; Double_t M=sqrt(L1*L2); // mutual inductance between solenoid and one coil/gamma . 0 < gamma < 1. // Double_t Rd = 0.06; // solenoid dump resistor, units are Ohms // old (incorrect) value, see next line. 9/21/2010 Double_t Rd = 0.072; // solenoid dump resistor, units are Ohms Double_t sigma = 0.0025; // half length of shorted region, units are m Double_t area = 4*sigma*sigma; // area of shorted region, units are m2 Double_t rho = 1.7e-8; // cu resistivity at room temperature, units are Ohm-m // Double_t rfactor = 2./75; // magneto-resistive effect x /1/RRR (cf Brindza e-mail of 4/30/2010 // Double_t RRR=100.; Double_t RRR=38.; // Double_t RRR=200.; Double_t rfactor = 1/RRR; Double_t Rs = rho*(2/75.)*2*sigma/area; // resistance of short, depends on resistivity at 4.4 K, length and area of shorted region. Double_t rf=10; Double_t Rs = Rs*rf; Double_t Rsolenoid = 1.1; // assume 1.1 m for radius of solenoid Double_t conductor_area=0.00762*0.00533; // assume 8 mm square conductor Double_t Rlcu = rho*rfactor*2*pi*Rsolenoid/conductor_area; //Dimension is Ohm*m , value taken from e-mail from Eugene. rho(2/75)6/0.008^2 is 4e-5 Ohm*m Double_t Rl_super = 1e-9; // use this value for "superconducting" to avoid numerical inversion problems. Double_t Rl = Rl_super; // superconducting Double_t I1_0 = 1500; //Value of current in solenoid loop at t=0, units are A Double_t I2_0 = 0; // Value of induced current in shorted loop at t=0, units are A printf ("L1=%g, L2=%g, M=%g, gamma0=%g, sigma=%g, area=%g, Rd=%g, Rs=%g, Rlcu = %g \n",L1,L2,M,gamma0,sigma,area,Rd,Rs,Rlcu); // Initial values for the currents in the solenoid and in the shorted loop Double_t I1[jmax]; Double_t I2[jmax]; Double_t I2abs[jmax]; Double_t I3[jmax]; Double_t Tshort[jmax]; Double_t Tloop[jmax]; Double_t t; Double_t tlimit=xmax; Double_t I1limit=0; Double_t I2limit=0; Double_t Tbath=4.2; Double_t Tshort[0]=Tbath; Double_t Tloop[0]=Tbath; Double_t kconst=2.5; // Thermal conductivity in W/cm-K Double_t ACmag=0.5; // Units are W/cm/K // Double_t ACmag=0.1; // Units are W/cm/K // Double_t ACmag=0.25; // Units are W/cm/K Double_t ACslope=0.09; // units are per deg K Double_t ACconst; // ACconst = ACmag*(1-ACslope*(T-Tbath)) superceeded by function; units are W/cm-K Double_t AreaCond=0.762*0.533; // Transverse area of conductor in units of cm2 x Dshort = volume to be cooled. Double_t ACwidth=0.533; // width of conductor Double_t Dshort=0.07; // length of short in cm (0.7 mm is thickness of ss band) Double_t Dloop=691; // length of loop in cm (radius=110 cm) Double_t ACPar[3]; Double_t ACPar[0]=ACmag; Double_t ACPar[1]=ACslope; Double_t ACPar[2]=1; Double_t A = L1; Double_t B = gamma*M; Double_t C = gamma*M; Double_t D = L2; Double_t temp=Tbath; // compute power dissipated Double_t Power1[jmax]; Double_t Power2[jmax]; Double_t Power3[jmax]; Double_t time[jmax]; Double_t Energy1=0; Double_t Energy2=0; Double_t Energy3=0; I1[0] = I1_0; I2[0] = I2_0; I2abs[0]=0; I3[0] = I1_0+I2_0; Power1[0]=Rd*I1_0*I1_0; Power2[0]=0; Power3[0]=0; // Loop over time interval, changing parameters at each step according to temperature. // Double_t deltat = (xmax-xmin)/jmax; Double_t deltat; Double_t logdeltat = (log(xmax)-log(xmin))/(jmax-1); // loop in logarithmic scale Double_t logt, logt0; printf ("Time deltat=%g\n",deltat); Double_t Tinput[1]; Double_t Parinput[4]={5.557161,-1.444021,1.281702,-0.213849}; Double_t Rpar[2]; Double_t Bfield=0; Double_t Rpar[0]=RRR; Double_t Rpar[1]=Bfield; for (j=1;j 0) { ACconst = transfer_func (Tinput, ACPar)*ACwidth/Tinput[0]; } else if { ACconst = 0; } // printf ("time=%f, Tinput=%f ACconst=%f\n",time,Tinput[0],ACconst); // determine kconst Tinput[0] = Tshort[j-1]; Double_t kconst = fit_heatcond (Tinput, Parinput); // printf ("Tinput[0]=%f, Parinput=%f %f %f %f, kconst=%f\n",Tinput[0],Parinput[0],Parinput[1],Parinput[2],Parinput[3],kconst); if (ACconst > 0) { Tshort[j] = Tbath + (APv/ACconst)*sqrt(ACconst/(AreaCond*kconst))*Dshort/2.; } else if { Tshort[j] = Tbath; } // printf ("time=%f APv=%f Tshort[j]=%f, ACconst=%f, kconst=%f \n",t,APv,Tshort[j],ACconst,kconst); // printf ("time=%f, Tshort=%f ACconst=%f\n",time,Tshort[j],ACconst); // use Eq. 18, GlueX-doc-1559 Double_t APv = (Power2[j-1]+Power3[j-1])/Dloop; Tinput[0] = Tloop[j-1] - Tbath + 0.01; ACPar[2] = 0; if (Tinput[0] > 0) { ACconst = transfer_func (Tinput, ACPar)*ACwidth/Tinput[0]; Tloop[j] = Tbath + APv/ACconst; } else if { ACconst = 0; Tloop[j] = Tbath; } // if (Tloop[j] < Tshort[j] && Rl > Rlcu/2) Tloop[j] = Tshort[j]; // printf ("time=%f, APv=%f, Tloop=%f ACconst=%f\n",t,APv,Tloop[j],ACconst); // change resistance based on current, temperature // printf ("j=%d, t=%f, I1[j-1]=%f, I2[j-1]=%f, I3[j-1]=%f\n",j,t,I1[j-1],I2[j-1],I3[j-1]); Tinput[0] = Tloop[j] ; // Tinput[0] = (Tloop[j] +Tshort[j])/2.; if (I3[j-1] > climit || Tshort[j] > 10) { // Rl = Rlcu; Rl = cu_resistance (Tinput,Rpar)*1e-8*2*pi*Rsolenoid/conductor_area; // reduce length of region that is assumed superconducting to 20 cm // Rl = cu_resistance (Tinput,Rpar)*1e-8*2*pi*Rsolenoid*0.1/conductor_area; // reduce length of region that is assumed superconducting to 0.1*2*pi*R cm // printf ("Note: Loosing Superconductivity, j=%d, time=%f, Tinput[0]=%f Tshort[j]=%f Tloop[j]=%f, RRR=%f Bfield=%f Rlcu=%f Rl=%f\n",j,t,Tinput[0],Tshort[j],Tloop[j],RRR,Bfield,Rlcu,Rl); } else if (I3[j-1] < climit && Tshort[j] < 7) { Rl = Rl_super; } else if (I3[j-1] > climit || Tloop[j] > 10) { Rl = cu_resistance (Tinput,Rpar)*1e-8*2*pi*Rsolenoid/conductor_area; // Rl = 1.72*rfactor*1e-8*2*pi*Rsolenoid/conductor_area; // assume constant resistance vs temperature // printf ("else j=%d, time=%f, Tinput[0]=%f Tshort[j]=%f Tloop[j]=%f, RRR=%f Bfield=%f Rlcu=%f Rl=%f\n",j,t,Tinput[0],Tshort[j],Tloop[j],RRR,Bfield,Rlcu,Rl); } // Compute the values of the resistances at each step in the loop based on current estimated temperature. // Rl and Rs depend on temperature Double_t B1 = Rd+Rl; Double_t B2= Rl; Double_t B3 = Rl; Double_t B4 = Rs+Rl; Double_t rlam1 = 0.5*(B1+B4) + 0.5*sqrt((B1-B4)*(B1-B4) + 4*B2*B3); Double_t rlam2 = 0.5*(B1+B4) - 0.5*sqrt((B1-B4)*(B1-B4) + 4*B2*B3); // printf ("B1=%g, B2=%g, B3=%g, B4=%g, rlam1=%g, rlam2=%g\n",B1,B2,B3,B4,rlam1,rlam2); Double_t Brlam1= B1 - rlam1; Double_t Brlam1s= Brlam1*Brlam1; Double_t Brlam2= B1 - rlam2; Double_t Brlam2s= Brlam2*Brlam2; Double_t Bs = B2*B2; // printf ("Brlam1=%g, Brlam2=%g, Brlam1s=%g, Brlam2s=%g, Bs=%g\n",Brlam1,Brlam2,Brlam1s,Brlam2s,B2); Double_t R1 = -B2/sqrt(Brlam1s + Bs); Double_t R2 = -B2/sqrt(Brlam2s + Bs); Double_t R3 = Brlam1/sqrt(Brlam1s+Bs); Double_t R4 = Brlam2/sqrt(Brlam2s+Bs); // printf (" R1=%g, R2=%g, R3=%g, R4=%g \n",R1,R2,R3,R4); // compute R inverse Double_t Rdet = R1*R4 - R3*R2; Double_t RI1 = R4/Rdet; Double_t RI2 = -R2/Rdet; Double_t RI3 = -R3/Rdet; Double_t RI4 = R1/Rdet; // printf (" RI1=%g, RI2=%g, RI3=%g, RI4=%g Rdet=%g\n",RI1,RI2,RI3,RI4,Rdet); Double_t RIM1 = RI1*A + RI2*C; Double_t RIM2 = RI1*B + RI2*D; Double_t RIM3 = RI3*A + RI4*C; Double_t RIM4 = RI3*B + RI4*D; // printf (" RIM: RIM1=%g, RIM2=%g, RIM3=%g, RIM4=%g \n",RIM1,RIM2,RIM3,RIM4); Double_t A1 = (RIM1*R1 + RIM2*R3)/rlam1; Double_t A2 = (RIM1*R2 + RIM2*R4)/rlam1; Double_t A3 = (RIM3*R1 + RIM4*R3)/rlam2; Double_t A4 = (RIM3*R2 + RIM4*R4)/rlam2; Double_t lam1 = 0.5*(A1+A4) + 0.5*sqrt((A1-A4)*(A1-A4) + 4*A2*A3); Double_t lam2 = 0.5*(A1+A4) - 0.5*sqrt((A1-A4)*(A1-A4) + 4*A2*A3); // printf("\n Use right eigenvalues\n"); // printf ("A1=%g, A2=%g, A3=%g, A4=%g, lam1=%g, lam2=%g\n",A1,A2,A3,A4,lam1,lam2); Double_t Alam1= A1 - lam1; Double_t Alam1s= Alam1*Alam1; Double_t Alam2= A1 - lam2; Double_t Alam2s= Alam2*Alam2; Double_t As = A2*A2; // printf ("Alam1=%g, Alam2=%g, Alam1s=%g, Alam2s=%g, As=%g\n",Alam1,Alam2,Alam1s,Alam2s,As); Double_t V1 = -A2/sqrt(Alam1s + As); Double_t V2 = -A2/sqrt(Alam2s + As); Double_t V3 = Alam1/sqrt(Alam1s+As); Double_t V4 = Alam2/sqrt(Alam2s+As); // printf (" V1=%g, V2=%g, V3=%g, V4=%g \n",V1,V2,V3,V4); // compute V inverse Double_t Vdet = V1*V4 - V3*V2; Double_t VI1 = V4/Vdet; Double_t VI2 = -V2/Vdet; Double_t VI3 = -V3/Vdet; Double_t VI4 = V1/Vdet; // printf (" VI1=%g, VI2=%g, VI3=%g, VI4=%g Vdet=%g\n",VI1,VI2,VI3,VI4,Vdet); Double_t VIM1 = VI1*A1 + VI2*A3; Double_t VIM2 = VI1*A2 + VI2*A4; Double_t VIM3 = VI3*A1 + VI4*A3; Double_t VIM4 = VI3*A2 + VI4*A4; // printf (" VIM: VIM1=%g, VIM2=%g, VIM3=%g, VIM4=%g \n",VIM1,VIM2,VIM3,VIM4); Double_t RV1 = R1*V1+R2*V3; Double_t RV2 = R1*V2+R2*V4; Double_t RV3 = R3*V1+R4*V3; Double_t RV4 = R3*V2+R4*V4; Double_t det = RV1*RV4 - RV3*RV2; Double_t RVI1 = VI1*RI1+VI2*RI3; Double_t RVI2 = VI1*RI2+VI2*RI4; Double_t RVI3 = VI3*RI1+VI4*RI3; Double_t RVI4 = VI3*RI2+VI4*RI4; // printf ("RV1=%g, RV2=%g, RV3=%g, RV4=%g\n",RV1,RV2,RV3,RV4); // printf ("RVI1=%g, RVI2=%g, RVI3=%g, RVI4=%g\n",RVI1,RVI2,RVI3,RVI4); Double_t Ip1_0 = (RV4*I1[j-1] - RV2*I2[j-1])/det; Double_t Ip2_0 = (-RV3*I1[j-1] + RV1*I2[j-1])/det; Double_t Ip1 = Ip1_0 * exp (-deltat/lam1); Double_t Ip2 = Ip2_0 * exp(-deltat/lam2); // printf ("j=%d,t=%g,I1[j-1]=%g,I2[j-1]=%g,Ip1=%g, Ip2=%g, det=%g\n",j,t,I1[j-1],I2[j-1],Ip1, Ip2,det); time[j] = t; I1[j] = RV1*Ip1 + RV2*Ip2; I2[j] = RV3*Ip1 + RV4*Ip2; I2abs[j] = abs(I2[j]); Power1[j] = I1[j]*I1[j]*Rd; Energy1 = Energy1 + Power1[j]*deltat; Power2[j] = I2[j]*I2[j]*Rs; Energy2 = Energy2 + Power2[j]*deltat; I3[j] = I1[j] + I2[j]; Power3[j] = I3[j]*I3[j]*Rl; Energy3 = Energy3 + Power3[j]*deltat; // printf("time=%g, current1=%g, current2=%g, current3=%g, Power1=%g, Power2=%g, Power3=%g\n",time[j],I1[j],I2[j],I3[j],Power1[j],Power2[j],Power3[j]); // printf("time=%g, Power1=%g, Power2=%g, Power3=%g\n",time[j],Power1[j],Power2[j],Power3[j]); // end loop over time go back for the next iteration. } // end of jloop Double_t Estored=0.5*L1*I1_0*I1_0; printf ("Estored=%g, Energy1=%g, Energy2=%g, Energy3=%g\n",Estored,Energy1,Energy2,Energy3); // TCanvas *c1 = new TCanvas("c1","c1 currents_cool",200,10,700,700); c1->SetGridx(); c1->SetGridy(); c1->SetBorderMode(0); c1->SetFillColor(0); c1->Divide(2,2); // c1->SetLogy(); TLegend *leg = new TLegend(0.60,0.75,0.95,0.9); c1->cd(1); c1_1->SetGridx(); c1_1->SetGridy(); c1_1->SetBorderMode(0); c1_1->SetFillColor(0); c1_1->SetLogx(); // c1_1->SetLogy(); Double_t ymin=-3000; Double_t ymax=10000; TGraph *dummy2 = new TGraph (npts,dummyx,dummyy); dummy2->SetMarkerColor(1); dummy2->SetMarkerStyle(21); dummy2->SetTitle(""); dummy2->GetXaxis()->SetRangeUser(xmin,xmax); dummy2->GetYaxis()->SetRangeUser(ymin,ymax); dummy2->GetXaxis()->SetTitleSize(0.04); dummy2->GetYaxis()->SetTitleSize(0.04); dummy2->GetYaxis()->SetTitleOffset(1.5); dummy2->GetXaxis()->SetTitle("Time (s)"); dummy2->GetYaxis()->SetTitle("Current (A)"); dummy2->GetXaxis()->SetNdivisions(505); dummy2->Draw("Ap"); sprintf (string,"Solenoid current\n"); TGraph *I1a = new TGraph (jmax,time,I1); leg->AddEntry(I1a,string,"p"); I1a->SetMarkerColor(1); I1a->SetMarkerStyle(21); I1a->SetMarkerSize(0.35); I1a->Draw("samep"); sprintf (string,"Short current\n",gamma); TGraph *I2a = new TGraph (jmax,time,I2); leg->AddEntry(I2a,string,"p"); I2a->SetMarkerColor(4); I2a->SetMarkerStyle(21); I2a->SetMarkerSize(0.35); I2a->Draw("samep"); TGraph *Isuma = new TGraph (jmax,time,I3); sprintf (string,"Loop current\n"); leg->AddEntry(Isuma,string,"p"); Isuma->SetMarkerColor(2); Isuma->SetMarkerStyle(21); Isuma->SetMarkerSize(0.35); Isuma->Draw("samep"); leg->Draw(); sprintf (string,"L_{1}=%g H\n",L1); t1 = new TLatex(0.20,0.84,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw(); sprintf (string,"L_{2}=%g H\n",L2); t1 = new TLatex(0.20,0.80,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw(); sprintf (string,"M=%g H\n",gamma*M); t1 = new TLatex(0.20,0.76,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw(); sprintf (string,"#gamma=%g\n",gamma); t1 = new TLatex(0.20,0.72,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw(); sprintf (string,"R_{dump}=%g #Omega\n",Rd); t1 = new TLatex(0.20,0.68,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw(); sprintf (string,"R_{short}=%g #Omega\n",Rs); t1 = new TLatex(0.20,0.64,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw(); sprintf (string,"R_{loop}=%g #Omega\n",Rlcu); t1 = new TLatex(0.20,0.60,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw(); sprintf (string,"#lambda_{1}=%g s\n",lam1); t1 = new TLatex(0.20,0.56,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw(); sprintf (string,"#lambda_{2}=%g s\n",lam2); t1 = new TLatex(0.20,0.52,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw(); sprintf (string,"Critical current=%.0f A\n",climit); t1 = new TLatex(0.20,0.48,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw(); TLegend *leg1 = new TLegend(0.60,0.75,0.95,0.9); c1->cd(2); c1_2->SetGridx(); c1_2->SetGridy(); c1_2->SetBorderMode(0); c1_2->SetFillColor(0); c1_2->SetLogx(); c1_2->SetLogy(); Double_t ymin=100; Double_t ymax=10000; TGraph *dummy2 = new TGraph (npts,dummyx,dummyy); dummy2->SetMarkerColor(1); dummy2->SetMarkerStyle(21); dummy2->SetTitle(""); dummy2->GetXaxis()->SetRangeUser(xmin,xmax); dummy2->GetYaxis()->SetRangeUser(ymin,ymax); dummy2->GetXaxis()->SetTitleSize(0.04); dummy2->GetYaxis()->SetTitleSize(0.04); dummy2->GetYaxis()->SetTitleOffset(1.5); dummy2->GetXaxis()->SetTitle("Time (s)"); dummy2->GetYaxis()->SetTitle("Current (A)"); dummy2->GetXaxis()->SetNdivisions(505); dummy2->Draw("Ap"); sprintf (string,"Solenoid current\n"); // TGraph *I1a = new TGraph (jmax,time,I1); leg1->AddEntry(I1a,string,"p"); I1a->SetMarkerColor(1); I1a->SetMarkerStyle(21); I1a->SetMarkerSize(0.35); I1a->Draw("samep"); sprintf (string,"Abs(Short current)\n",gamma); TGraph *I2b = new TGraph (jmax,time,I2abs); leg1->AddEntry(I2b,string,"p"); I2b->SetMarkerColor(4); I2b->SetMarkerStyle(21); I2b->SetMarkerSize(0.35); I2b->Draw("samep"); TGraph *Isuma = new TGraph (jmax,time,I3); sprintf (string,"Loop current\n"); leg1->AddEntry(Isuma,string,"p"); Isuma->SetMarkerColor(2); Isuma->SetMarkerStyle(21); Isuma->SetMarkerSize(0.35); Isuma->Draw("samep"); leg1->Draw(); c1->cd(3); c1_3->SetGridx(); c1_3->SetGridy(); c1_3->SetBorderMode(0); c1_3->SetFillColor(0); c1_3->SetLogx(); c1_3->SetLogy(); // Double_t xmin=0.1; // Double_t xmax=1000; Double_t ymin=0.01; Double_t ymax=200000; TGraph *dummy = new TGraph (npts,dummyx,dummyy); dummy->SetMarkerColor(1); dummy->SetTitle(""); dummy->GetXaxis()->SetRangeUser(xmin,xmax); dummy->GetYaxis()->SetRangeUser(ymin,ymax); dummy->GetXaxis()->SetTitleSize(0.04); dummy->GetYaxis()->SetTitleSize(0.04); dummy->GetYaxis()->SetTitleOffset(1.5); dummy->GetXaxis()->SetTitle("Time (s)"); dummy->GetYaxis()->SetTitle("Power (W)"); dummy->GetXaxis()->SetNdivisions(505); dummy->Draw("Ap"); TGraph *P1 = new TGraph (jmax,time,Power1); TLegend *leg3 = new TLegend(0.60,0.75,0.95,0.9); P1->SetMarkerColor(1); P1->SetMarkerStyle(21); P1->SetMarkerSize(0.35); P1->Draw("samep"); sprintf (string,"Stored Energy=%g J\n",Estored); t1 = new TLatex(0.20,0.80,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.04); t1->Draw(); sprintf (string,"Energy Solenoid=%g J\n",Energy1); t1 = new TLatex(0.20,0.75,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.04); t1->Draw(); sprintf (string,"Energy Short=%g J\n",Energy2); t1 = new TLatex(0.20,0.70,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.04); t1->Draw(); sprintf (string,"Energy Loop=%g J\n",Energy3); t1 = new TLatex(0.20,0.65,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.04); t1->Draw(); TGraph *P2 = new TGraph (jmax,time,Power2); P2->SetMarkerColor(4); P2->SetMarkerStyle(21); P2->SetMarkerSize(0.35); P2->Draw("Samep"); TGraph *P3 = new TGraph (jmax,time,Power3); P3->SetMarkerColor(2); P3->SetMarkerStyle(21); P3->SetMarkerSize(0.35); P3->Draw("Samep"); c1->cd(4); c1_4->SetGridx(); c1_4->SetGridy(); c1_4->SetBorderMode(0); c1_4->SetFillColor(0); c1_4->SetLogx(); c1_4->SetLogy(); // Double_t xmin=0.1; // Double_t xmax=1000; Double_t ymin=1; Double_t ymax=1000; TGraph *dummy = new TGraph (npts,dummyx,dummyy); dummy->SetMarkerColor(1); dummy->SetTitle(""); dummy->GetXaxis()->SetRangeUser(xmin,xmax); dummy->GetYaxis()->SetRangeUser(ymin,ymax); dummy->GetXaxis()->SetTitleSize(0.04); dummy->GetYaxis()->SetTitleSize(0.04); dummy->GetYaxis()->SetTitleOffset(1.5); dummy->GetXaxis()->SetTitle("Time (s)"); dummy->GetYaxis()->SetTitle("Temperature (K)"); dummy->GetXaxis()->SetNdivisions(505); dummy->Draw("Ap"); sprintf (string,"ACmag=%.3f W/cm/K\n",ACmag); t1 = new TLatex(0.20,0.84,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw(); sprintf (string,"ACslope=%.3f K^{-1}\n",ACslope); t1 = new TLatex(0.20,0.80,string); t1->SetTextColor(1); t1->SetNDC(); t1->SetTextSize(0.03); t1->Draw(); TGraph *T2 = new TGraph (jmax,time,Tshort); TLegend *leg4 = new TLegend(0.60,0.75,0.95,0.9); sprintf (string,"Short Temperature\n"); leg4->AddEntry(T2,string,"p"); T2->SetMarkerColor(4); T2->SetMarkerStyle(21); T2->SetMarkerSize(0.35); T2->Draw("samep"); TGraph *T3 = new TGraph (jmax,time,Tloop); sprintf (string,"Loop Temperature\n"); leg4->AddEntry(T3,string,"p"); T3->SetMarkerColor(2); T3->SetMarkerStyle(21); T3->SetMarkerSize(0.35); T3->Draw("samep"); leg4->Draw(); // save canvas to file Int_t ig = gamma*1000; Int_t ir = rf*10; sprintf(filename,"currents_cool_c1_%d_%d.eps",ig,ir); c1->SaveAs(filename); sprintf(filename,"currents_cool_c1_%d_%d.png",ig,ir); c1->SaveAs(filename); } Double_t fit_heatcond (Double_t *x, Double_t *par) { // Fit ln(thermal conductivity) as a funciton of ln(Temperature). Returns value in W/cm-K Double_t p0=par[0]; Double_t p1=par[1]; Double_t p2=par[2]; Double_t p3=par[3]; Double_t x1=x[0]; Double_t pi=3.14159; char string[256]; Double_t func; if (x1 < 4) { x1 = 4; } elseif (x1 > 70) { x1 = 70; } x1 = log(x1); func = p0 + p1*x1 + p2*x1*x1 + p3*x1*x1*x1; /*sprintf (string,"x1=%f func=%f\n",x1,func); printf ("string=%s",string);*/ if (func < 50) { return exp(func)/100.; // change to W/cm-K } else { return exp(50)/100.; // change to W/cm-K } } Double_t transfer_func (Double_t *x, Double_t *par) { // parameterization of cooling over a large temperature range due to nucleate and film boiling. // note that ACconst = heat transfer q (W/cm2)*deltax/(Delta T) Double_t ACmag=par[0]; Double_t ACslope=par[1]; Double_t option=par[2]; // option=1 upper curve, option=0 lower curve Double_t x1=x[0]; // Note that input temperature should be given relative to the helium bath Double_t pi=3.14159; Double_t Tbath=4.2; Double_t T1=1; // degrees relative to bath temperature Double_t T2=10; // degrees relative to bath temperature Double_t deltax=0.533; // width of conductor (cm) char string[256]; if (x1 < 0.001) { return 0; } Double_t ACconst; Double_t factor; if (option) { factor = 1; } else if { factor=10; } if (x1 < T1/factor) { ACconst = ACmag*(1-ACslope*T1); } else if (x1 < T2/factor) { ACconst = ACmag*(1-ACslope*T1)/(x1*factor); } else { ACconst = ACmag*(1-ACslope*T2); } Double_t q = ACconst*x1/deltax; // note that ACconst * Delta T = heat transfer q (W/cm)/deltax = q(W/cm2), p. 202 of Helium Cryogenics by Van Sciver /*sprintf (string,"x1=%f func=%f\n",x1,func); printf ("string=%s",string);*/ return q; } Double_t cu_magnetoresistivity (Double_t *x, Double_t *par) { // Parameterization of the magneto resistance of Cu from Benz Journal of Applied Physics vol 40, number 1 April 1969, p. 2003 Double_t S0=par[0]; // Resistivity ratio at H=0; // Double_t x1=x[0]*1e-6; // independent variable is H*S0 x 10^-6 G Eq. (3). Double_t x1=x[0]*S0*1e-2; // independent variable is H*S0 x 10^-6 G Eq. (3), Use H in Tesla (10^4 G) Double_t pi=3.14159; Double_t R0=1.72/S0; // Units are 10^-8 Ohm-m; Resistivity at 4.2 deg K. Double_t k1=2.368e-4; // Constants in polynomial fit to data X<2 Double_t k2=7.389e-2; Double_t k3=6.505e-2; Double_t k4= -7.121e-3; Double_t c1=-0.1649; // Constants in polynomial fit to data X>2 Double_t c2=0.2603; Double_t c3=-1.129e-3; Double_t c4= 4.011e-6; char string[256]; Double_t func=0; if (x1 < 2) { func = k1 + k2*x1 + k3*x1*x1 + k4*x1*x1*x1; } else if (x1 < 130) { func = c1 + c2*x1 + c3*x1*x1 + c4*x1*x1*x1; } // Return value of resistivity instead of fractional change from H=0; // func = R0*(1+func); // units are 10^-8 Ohm-m; /*sprintf (string,"x[0]=%f, x1=%f S0=%f, R0=%f, func=%f\n",x[0],x1,S0, R0,func); printf ("string=%s",string);*/ return func; } Double_t cu_resistance (Double_t *x, Double_t *par) { // x[0] is the temperature, cu_resistance is in 10^-8 Ohm-m. // par[0] is the RRR // par[1] is the B field in Tesla // Fit ln(electrical resistivity) as a funciton of ln(Temperature) for copper // Use Fermi function to fit resistivity vs temperature // Note that function fits data best for RRR ~ 20 - 200, as the scaled function is scaled from a RRR=42, but the data dependence // change slightly as a function of RRR. Double_t rho=4.202; Double_t slope=1.61; Double_t cc2=-0.6352; Double_t cc3=0.4155; Double_t scale=1; Double_t mean=log(75); // Double_t RRR=par[0]; Double_t RRR=par[0]*1.72/1.545; //Note that RRR should be referenced to T=273K. (not 300 K). Adjust for difference Double_t B[1]; B[0] = par[1]; // pass on the value of the B-field. // scale RRR by the B field factor. Double_t deltaB = cu_magnetoresistivity (B, par); RRR = RRR/(1 + deltaB); Double_t Rcu_300=1.72; // in units of 10^-8 Ohm-m Double_t x1=log(x[0]); Double_t pi=3.14159; /*Double_t RRR850= 850; // use this RRR value as reference. Double_t logRRR850 = log(Rcu_300/RRR850);*/ Double_t RRR42= 42; // use this RRR value as reference. Double_t logRRR42 = log(Rcu_300/RRR42); char string[256]; Double_t func; Double_t arg; arg = -(slope*(x1-mean)+cc2*(x1-mean)*(x1-mean)+cc3*(x1-mean)*(x1-mean)*(x1-mean) ); if (arg < -100) { func = rho; } else if (arg > 100) { func = logRRR42; } else { func = logRRR42 + rho/(1+exp(arg)); } // reference to room temperature scale = log(RRR)/log(RRR42); // printf ("x[0]=%f RRR %f RRR42 %f scale=%f \n",x[0],RRR,RRR42,scale); Double_t scaled = exp(log(Rcu_300)-scale*(log(Rcu_300)-func)); // printf("x[0] %f B[0] %f par[0] %f delta B %f\n",x[0],B[0],par[0],deltaB); return scaled; }