c c exportrates.f c c This code calculates neutrino capture rates and the uncertainties c in neutrino fluxes and capture rates. c c Revised 10/08/2004. c c *************************** c Precepts and numerical data c *************************** c c The precepts and values embodied in this subroutine are described in a c series of RMP papers: J. N. Bahcall, RMP 50, 881 (1978); Bahcall et al. c RMP, 54, 767 (1982)--see especially Eq. (23), pg. 787; Bahall and R. K. c Ulrich, RMP 60, 297 (1988)--see especially Table XV and Eq. (16), pg. c 330 and Eq. (17), pg 333; Bahcall and M. P. Pinsonneault, RMP 64, 885 c (1992) and RMP 67, 781 (1995). There are some more recent neutrino c cross sections for gallium (J. N. Bahcall, Phys Rev C, 56, 3391, 1997), c and for chlorine (J. N. Bahcall et al., Phys Rev C, 54, 411, 1996). c The Li absorption cross section for 7Be neutrinos is from Bahcall, Phys c Rev D, 49, 3923 (1994). The earlier neutrino cross sections are from c Bahcall ``Neutrino Astrophysics'', 1989, unless specified otherwise. c Two recent discussions of rates and their uncertainties are Bahcall, c S. Basu, and M. H. Pinsonneault, Physics Letters B, 433, 1-8 (1998), c see especially Table 2, and J. N. Bahcall, M. H. Pinsonneault, Astroph. c Journal 555, 990-1012 (July 10, 2001), astro-ph/0010346, see especially c Section 2 and Table 7. The most recent discussion of the nuclear rates is c in Bahcall and Pinsonneault, Physical Review Letters, 92, 121301 (2004), c astro-ph/0402114. c c Input uncertainties are 3-sigma; output uncertainties are 1-sigma. (Gallium c neutrino absorption cross sections are an exception.) c c Exportrates.dat contains the input neutrino fluxes and the uncertainties. c c Can read in several different cases with different specified uncertainties. c This allows one to take account of assymetric errors. c c c c ***************************************** c Information re: order of neutrino sources c ***************************************** c c The order of the neutrino sources is: pp,pep,Be7,B8, N13,O15,F17,hep. c c************ c Definitions c************ c c Stndrd(I): Standard solar neutrino fluxes. c c c CrossCl(I): Cross sections for the chlorine detector. Unit: 10**-46 c cm**2. c c CrossGa(I): Cross sections for the gallium detector. Unit: 10**-46 c cm**2. c c CrossLi(I): Cross sections for the Li detector. Unit: 10**-46 c cm**2. c c UncrossCl(I) Fractional uncertainty in Cl neutrino absorption cross c section. Note cross sectional uncertainties are fractional c uncertainties. c c c UncrossGa(I) Fractional uncertainty in Ga neutrino absorption cross c section. c c UncrossLi(I) Fractional uncertainty in Li neutrino absorption cross c section. c c c UnlowerGa(I) Same as UncrossGa(I) except to be used for calculating c lower limits of uncertainties. c c c PartCl(I): Product of flux times cross section for Ith neutrino c source. c c PartGa(I): Same as PartCl(I) except for Ga instead of Cl. c c c PartLi(I): Product of flux times cross section for Ith neutrino c source. c c UncertJ: The uncertainty in the parameter J defined by c [1 + (dx sub j)/(x sub j) ] as in equation 23 c Bahcall-Ulrich RMP '82. If the 3sigma error in c x is 3*epsilon*x, then UncertaintJ = 1 =3*epsilon, c dx = 3*epsilon*x . c c Uncertsq(J): Squares of uncertainties: Uncert(J)*Uncert(J). c c opc(I): Opacity uncertainty. The total fractional differences c in Table 7 of Bahcall and Pinsonnault (1992). c The differences are computed using the first and third c rows of Table 7, comparing the Los Alamos and the much c improved Livermore opacities. c Note: DelopCl = (Sum over I)[PartCl(I)*opc(I)] c c diffuse(I): [Flux(I)_(without diffusion) - Flux(with diffusion)]/ c Flux(with diffusion). Values taken from Models 9 and 10 c of Bahcall-Pinsonneault 1995. c c S11(I): Partial Derivatives of neutrino fluxes with respect to p-p c rate. From Bahcall and Ulrich 88 c c S34(I): Partial Derivatives of neutrino fluxes with respect to 3-4 c rate. From Bahcall and Ulrich 88. c c S33(I): Partial Derivatives of neutrino fluxes with respect to 3-3 c rate. From Bahcall and Ulrich 88. c c Shep(I): Partial Derivatives of neutrino fluxes with respect to hep c cross section. Obvious because hep is so rare a reaction. c c S114(I): Partial Derivatives of neutrino fluxes with respect to c N14-p rate. From Bahcall and Ulrich 88. c c Sunlum(I): Partial Derivatives of neutrino fluxes with respect to the c solar luminosity. From Bahcall and Ulrich 88. c c ZoverX(I): Partial derivatives of neutrino fluxes with respect to the c ratio of Z/X . From Bahcall and Ulrich 88 or BP04. Can choose. c c Age(I) Partial derivatives of neutrino fluxes with respect to the c solar age. From Bahcall and Ulrich 88. c c UncertComposition(I) The uncertainties in the neutrino fluxes and rates due to c individually computed heavy element abundances. The first 8 values c are uncertainties of solar neutrino fluxes. The next three values are c the uncertainties of the Cl, Ga, and Li event rates, respectively. All c of these uncertainties are read in as 1 sigma uncertainties. They were c calculated with the aid of partial derivatives of each flux with respect c to each of the principal heavy element abundances. c c IZ: Use individual composition uncertainties if IZ = 7. Otherwise, uses c historical average for composition uncertainties. The individual c composition uncertainties are obtained from compositionuncertainties.f c and are 1 sigma uncertainties. c c*************** c Uncertainties c*************** c c Nuclear physics Uncertainties are from the RMP Solar Fusion paper c Adelberger et al. 98, Rev. Mod. Phys. 70, 1265, unless explicitly c noted as updated. c c Z/X: Use the uncertainty quoted by Bahcall and Pinsonneault (1995), c Eq. (7). This is the 'historical' option and has been replaced in c this version by the individually calculated uncertainties due to each c heavy element abundance. c c The 17F flux is proportional to the 16O(p,\gamma)17F reaction rate. This c uncertainty is added in quadrature to the other 17F flux uncertainties. c Since this reaction affects only the 17F flux, its effect is only c included in the 17F uncertainty calculation. The Adelberger et al. c uncertainty is adopted. c c Age uncertainty is based upon Wasserburg analysis (see appendix of BP 95). c c Uncertdiffuse is the fractional uncertainty in the rate of diffusion. c The total effect of diffusion is determined by comparing Models 9 and 10 c of Table 3 of Bahcall-Pinsonneault 1995. c c UncertBe7e is the fractional uncertainty (3 sigma) in the Be7 + e c capture rate, determined by the 1997 paper of Gruzinov and Bahcall. c Ap J 490, 437. This uncertainty only affects (inversely) the 8B c neutrino flux since the 7Be electron capture is so much faster c than 7Be proton capture. c c UnGalower is the 1 sigma lower limit for the gallium rate, computed with c lower-limit neutrino absorption cross section uncertainties for Ga. c OneSigmaGa is the 1 sigma upper limit for the gallium rate, computer with c upper-limit neutrino absorption cross section uncertainties for Ga. c LowerGanucr is the 1 sigma uncertainty in the Ga rate due just to neutrino c cross sections. In the code, it is referred to as DellowerGa and can be c printed out with Write(1,122) statement if desired. c c********** c Averages: c********** c c The variables av**** are used to average uncertainties across c cases. Each variable is indexed as defined by the parameters c above. Note that this means some locations are not used. (avnucr c is only used for Cl, Ga and Li for example.) However, this results c in simpler code and does not waste much array space. c c********************* c Dimension Statements c********************* c Dimension Stndrd(8) Dimension UncrossCl(8),UncrossGa(8),UncrossLi(8),UnlowerGa(8) Dimension CrossCl(8),CrossGa(8),CrossLi(8) Dimension S11(8),S34(8),S114(8),S17(8) Dimension S33(8),Shep(8),Sunlum(8),Age(8) Dimension PartCl(8),PartGa(8),PartLi(8) Dimension opc(8) Dimension ZoverX(8),diffuse(8) Dimension Uncert(8), Uncertsq(8) Dimension av11(11),av33(11),av34(11),av114(11),av17(11),avhep(11) Dimension avZX(11), avopac(11), avlum(11), avage(11) Dimension avdiffuse(11),avBe7e(11),avnucr(11),avUncert(11) Dimension UncertComposition(11) c c****************** c Type declarations ******************* c Real N1311,N1333,N1334,N1317,N13ZX,N13opac,N13114,N13age,N13lum, * N13diffuse Real O1511,O1533,O1534,O1517,O15ZX,O15opac,O15114,O15age,O15lum Integer PP,PEP,BE7,B8,N13,O15,F17,HEP,CL,GA,LI Character*4 Source(11) c c ******************** c Parameter statements c ******************** c c Here are the parameters used for indexing. c Parameter (PP=1,PEP=2,BE7=3,B8=4,N13=5,O15=6) Parameter (F17=7,HEP=8,CL=9,GA=10,LI=11) c c***************************************************************** c Input neutrino absorption cross sections and their uncertainties. c***************************************************************** c c The cross section uncertainties given here are valid only IF nothing c changes the shape of the neutrino energy spectrum. Uncertainties as a c function of energy should be used if the neutrino energy spectrum is c changed. See Phys Rev C 56, 3391 (1997). c c The input uncertainties are 3 sigma fractional uncertainties for c Cl and Li and 1 sigma for Ga. c c Take 37Cl cross section for 8B neutrinos from Bahcall c et al. (1996), Phys. Rev. C,54, 411. c Data CrossCl/0.0,16.,2.4,11400.,1.7,6.8,6.9,42600./ Data UncrossCl/0.00,0.06,0.06,0.096,0.06,0.06,0.06,0.11/ c c Gallium cross sections are from Phys Rev C 56, 3391,1997, see c Table VI. These are 1-sigma uncertainties. c Data CrossGa/11.72,204.,71.7,24000.,60.4,113.7,113.9,71400./ Data UncrossGa/0.023,0.17,0.07,0.32,0.06,0.12,0.12,0.32/ Data UnlowerGa/0.023,0.07,0.03,0.15,0.03,0.06,0.06,0.16/ c c Convert gallium uncertainties to 3 sigma. c Do I = 1,8 UncrossGa(I) = 3.0*UncrossGa(I) UnlowerGa(I) = 3.0*UnlowerGa(I) End do c c Have included the Li cross section for 7Be from PRD,49,3923(1994). c Data CrossLi/0.0,655.,19.,39000.,42.4,246.,249.,84000./ Data UncrossLi/0.0,0.06,0.05,0.17,0.055,0.06,0.06,0.21/ c c****************************************************************** c Partial derivatives of neutrino fluxes with respect to parameters c****************************************************************** c c The tabulated quantites are the logarithmic partial derivatives c of neutrino fluxes with respect to input parameters, defined by c Eq. (7.3) of Neutrino Astrophysics. Thus we tabulate c parameter_{i} = \partial (ln (flux_i) )/\partial (ln parameter) . c c The values given here for opc(I) are taken from Table 7 of c Bahcall-Pinsonneault (1992). They are twice the difference c between rows 1 and 3, divided by the sum of rows 1 and 3 (entry c by entry). The defining equation is Eq. (20) of Bahcall-Pinsonneault c 1992. The other values are taken from Bahcall-Ulrich 1988, Table XV. c The values of diffuse are the fractional differences between c Models 9 and 10 of Table 3 of Bahcall-Pinsonneault 1995. The numbers c shown are [phi(NoDiff) - phi(Diff)]/phi(Diff). To get the 3-sigma c uncertainty due to diffusion, must multiply by Uncertdiffuse, which is c read in from exportrates.dat . c Data S11/0.14,-0.17,-.97,-2.59,-2.53,-2.93,-2.94,-0.08/ Data S33/0.03,0.05,-0.43,-0.40,0.02,0.02,0.02,-0.45/ Data S34/-0.06,-0.09,0.86,0.81,-0.05,-0.05,-0.05,-0.08/ Data S17/0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0/ Data S114/-0.02,-0.02,0.0,0.01,0.85,1.0,.01,-0.01/ Data Shep/0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0/ Data Sunlum/0.73,0.87,3.4,6.76,5.16,5.94,6.25,0.12/ c Data ZoverX/-0.08,-0.17,0.58,1.265,1.86,2.03,2.09,-0.22/ c c The following ZoverX values were computed for BP04. c Data ZoverX/-0.084,-0.171,0.619,1.364,1.897,2.056,2.169,-0.242/ c Data opc/0.0083,0.0142,-.0829,-.157,-.100,-.123, *-0.128,0.032/ Data age/-0.07,0.0,0.69,1.28,1.01,1.27,1.29,-0.11/ Data diffuse/0.01692,0.02857,-0.12039,-0.26737,-0.34142, *-0.36697,-0.37963,0.04959/ c c Here are the labels used for output. c Data Source/'pp','pep','7Be','8B','13N','15O','17F','hep','Cl', * 'Ga','Li'/ c c Set all av variables to 0 for summing for averages. c Data av11/11*0.0/, av33/11*0.0/, av34/11*0.0/, av114/11*0.0/ * av17/11*0.0/, avZX/11*0.0/, avopac/11*0.0/, avlum/11*0.0/, * avage/11*0.0/, avdiffuse/11*0.0/, avBe7e/11*0.0/, * avnucr/11*0.0/, avUncert/11*0.0/, avhep/11*0.0/ c c********************************* c Open the output and input files. c********************************* c Open(1,file = 'exportrates.output', status = 'unknown') c c Input neutrino fluxes. Units of 10^10. c Open(3,file = 'exportrates.dat', status = 'old') c c*********************************************** c Read in the neutrino fluxes and uncertainties. c*********************************************** c c First eight lines of exportrates.dat are not read; they are skipped. c Then read in the 8 neutrino fluxes, Stndrd(I) from line 9. c If IZ = 7,on row 10, then will use individual composition uncertainties. c The flux uncertainties are read from row 11 and are used only if IZ = 7. c Otherwise, proceed with different cases of uncertainties, using historical c average for composition uncertainties as expressed in terms of Z/X. Thus c IZ not equal to 7 reduces to method used to calculate abundance uncrtainties c prior to October 2004. c Read(3,300)Skip Read(3,*)(Stndrd(I), I = 1,8) Write(6,180) Write(1,180) Write(6,182) Write(1,182) Write(6,*)(Stndrd(I), I = 1,8) Write(1,*)(Stndrd(I), I = 1,8) Read(3,*)IZ Write(6,185)IZ Write(1,185)IZ Read(3,*)(UncertComposition(I),I = 1,11) c c Convert composition uncertainties to 3 sigma. Use only if IZ = 7. c23456789012345678901234567890123456789012345678901234567890123456789012 c c Write out individual composition uncertainties (1 sigma). c Write(6,181) Write(1,181) Write(6,184)(UncertComposition(I),I = 1,11) Write(1,184)(UncertComposition(I),I = 1,11) c Do I = 1,11 UncertComposition(I) = 3.0*UncertComposition(I) EndDo c c Write out individual composition uncertainties (3 sigma). c c Write(6,*) (UncertComposition(I), I = 1,11) c Write(1,*) (UncertComposition(I), I = 1,11) c c c********************************************************************** c Begin Do-loop over a number of different cases of the uncertainties. c Stop when the first uncertainty in a row of uncertainties is negative. c********************************************************************** c Ncase = 0 Do IUncert = 1, 10 Read(3,*)Uncert11,Uncert33,Uncert34,Uncert17,Uncert114, *UncertZX,Uncertlum,Uncertage,Uncertdiffuse,UncertBe7e,Uncerthep If(Uncert11 .lt. 0.000) go to 999 Ncase = Ncase + 1 c c Write the label for the top of the page and the case number, c IUncert (indicates which set of uncertainties are being used). c c Note: c Uncomment all write(1,400) statements to get top of page functionality c c Write(1,400) Write(1,100) Write(1,102) IUncert c c********************************************* c Calculate Rates: Fluxes times cross sections c********************************************* c RateCl = 0.0 RateGa = 0.0 RateLi = 0.0 c c Calculate each source's contribution to each experiment. c Do I = 1,8 PartCl(I) = Stndrd(I)*CrossCl(I) PartGa(I) = Stndrd(I)*CrossGa(I) PartLi(I) = Stndrd(I)*CrossLi(I) RateCl = RateCl + PartCl(I) RateGa = RateGa + PartGa(I) RateLi = RateLi + PartLi(I) End do c c Write out individual contributions to each target, but only for case 1. c If(IUncert.eq.1) then Write(1,103) Write(1,104)(Source(I),Stndrd(I), * PartCl(I),PartGa(I),PartLi(I),I=1,8) Endif c ************************************ c Begin calculating the uncertainties c ************************************ c c The equation used here is equation (23) of RMP 82. I calculate here c for j = p-p, 3-3, 3-4, 1-14,1-7, Z/X, luminosity, opacity, age, and c neutrino cross sections. I first calculate all of the partial c uncertainties (the sum over i for each parameter separately c in equation 23 ) and then sum them quadratically and take the square c root. At this stage, the uncertainties are still 3 sigma. c Del11Cl = 0.0 Del11Ga = 0.0 Del11Li = 0.0 Do I = 1,8 Del11Cl = Del11Cl + PartCl(I)*( Uncert11**S11(I) -1.0 ) Del11Ga = Del11Ga + PartGa(I)*( Uncert11**S11(I) -1.0 ) Del11Li = Del11Li + PartLi(I)*( Uncert11**S11(I) -1.0 ) End Do c Del33Cl = 0.0 Del33Ga = 0.0 Del33Li = 0.0 Do I = 1,8 Del33Cl = Del33Cl + PartCl(I)*( Uncert33**S33(I) -1.0 ) Del33Ga = Del33Ga + PartGa(I)*( Uncert33**S33(I) -1.0 ) Del33Li = Del33Li + PartLi(I)*( Uncert33**S33(I) -1.0 ) End Do c Del34Cl = 0.0 Del34Ga = 0.0 Del34Li = 0.0 Do I = 1,8 Del34Cl = Del34Cl + PartCl(I)*( Uncert34**S34(I) -1.0 ) Del34Ga = Del34Ga + PartGa(I)*( Uncert34**S34(I) -1.0 ) Del34Li = Del34Li + PartLi(I)*( Uncert34**S34(I) -1.0 ) End Do c Del114Cl = 0.0 Del114Ga = 0.0 Del114Li = 0.0 Do I = 1,8 Del114Cl = Del114Cl + PartCl(I)*( Uncert114**S114(I) -1.0 ) Del114Ga = Del114Ga + PartGa(I)*( Uncert114**S114(I) -1.0 ) Del114Li = Del114Li + PartLi(I)*( Uncert114**S114(I) -1.0 ) End Do c c Only the B8 neutrino flux is affected by an uncertainty in S_{17} . c (The formula for the uncertainty reduces to this simple case if the c partial derivatives are all either zero except for the one which is c equal to unity.) c Del17Cl = 0.0 Del17Ga = 0.0 Del17Li = 0.0 Del17Cl = PartCl(4)*( Uncert17 -1.0 ) Del17Ga = PartGa(4)*( Uncert17 -1.0 ) Del17Li = PartLi(4)*( Uncert17 -1.0 ) c c*************************** c Uncertainty due to S_(hep) c*************************** c c Only the hep neutrino flux is affected by an uncertainty in S_{hep} . c DelhepCl = 0.0 DelhepGa = 0.0 DelhepLi = 0.0 DelhepCl = PartCl(8)*( Uncerthep -1.0 ) DelhepGa = PartGa(8)*( Uncerthep -1.0 ) DelhepLi = PartLi(8)*( Uncerthep -1.0 ) c c******************************************************* c Uncertainty due to heavy element composition of the Sun c******************************************************** c DelZXCl = 0.0 DelZXGa = 0.0 DelZXLi = 0.0 c If (IZ .eq. 7) Then c C Individual compositions c DelZXCl = UncertComposition(9) DelZXGa = UncertComposition(10) DelZXLi = UncertComposition(11) Else c C Z/X uncertainty from historical average c Do I =1,8 DelZXCl = DelZXCl + PartCl(I)*( UncertZX**ZoverX(I) -1.0 ) DelZXGa = DelZXGa + PartGa(I)*( UncertZX**ZoverX(I) -1.0 ) DelZXLi = DelZXLi + PartLi(I)*( UncertZX**ZoverX(I) -1.0 ) End Do Endif c c***************************************** c Uncertainty due to luminosity of the Sun c***************************************** c DellumCl = 0.0 DellumGa = 0.0 DellumLi = 0.0 Do I = 1,8 DellumCl = DellumCl + PartCl(I)*( Uncertlum**Sunlum(I) -1.0 ) DellumGa = DellumGa + PartGa(I)*( Uncertlum**Sunlum(I) -1.0 ) DellumLi = DellumLi + PartLi(I)*( Uncertlum**Sunlum(I) -1.0 ) End Do c DelageCl = 0.0 DelageGa = 0.0 DelageLi = 0.0 Do I = 1,8 DelAgeCl = DelAgeCl + PartCl(I)*( UncertAge**Age(I) -1.0 ) DelAgeGa = DelAgeGa + PartGa(I)*( UncertAge**Age(I) -1.0 ) DelAgeLi = DelAgeLi + PartLi(I)*( UncertAge**Age(I) -1.0 ) End Do c c************************************* c Calculate uncertainty due to opacity c************************************* c DelopcCl = 0.0 DelopcGa = 0.0 DelopcLi = 0.0 Do I = 1,8 DelopcCl = DelopcCl + PartCl(I)*opc(I) DelopcGa = DelopcGa + PartGa(I)*opc(I) DelopcLi = DelopcLi + PartLi(I)*opc(I) End Do c c********************************************************************** c Calculate the diffusion uncertainty assuming 3-sigma = no diffusion. c********************************************************************** c c Uncertdiffuse is read in from exportrates.dat. c DeldiffuseCl = 0.0 DeldiffuseGa = 0.0 DeldiffuseLi = 0.0 Do I = 1,8 DeldiffuseCl = DeldiffuseCl + PartCl(I)*diffuse(I)*Uncertdiffuse DeldiffuseGa = DeldiffuseGa + PartGa(I)*diffuse(I)*Uncertdiffuse DeldiffuseLi = DeldiffuseLi + PartLi(I)*diffuse(I)*Uncertdiffuse End Do c c*************************************************** c Uncertainty in neutrino absorption cross sections. c*************************************************** c c For gallium only (so far), I have explicitly inclucded of the calculation c of the lower limit of neutrino absorption cross section uncertainties. c This is represented by UnlowerGa(I). For Cl and Li, I have so far c used only symmetric uncertainties. c c Will keep the high energy neutrinos separate from the lower energy c neutrinos. It is reasonable from a nuclear physics point of view c that the uncertainties in the two separate groups (high and low energy) c should be coherent, but that they should be added separately. c c Will now break up into high and low energy parts. c The high energy neutrinos go into group 1 and the lower energy c (all but B8 and hep) go into group 2. c DelnucrCl1 = 0.0 DelnucrGa1 = 0.0 DellowerGa1 = 0.0 DelnucrLi1 = 0.0 DelnucrCl2 = 0.0 DelnucrGa2 = 0.0 DellowerGa2 = 0.0 DelnucrLi2 = 0.0 Do I = 1,8 If(I.eq.4 .or. I.eq.8) Then DelnucrCl1 = DelnucrCl1 + PartCl(I)*UncrossCl(I) DelnucrGa1 = DelnucrGa1 + PartGa(I)*UncrossGa(I) DellowerGa1 = DellowerGa1 + PartGa(I)*UnlowerGa(I) DelnucrLi1 = DelnucrLi1 + PartLi(I)*UncrossLi(I) Else DelnucrCl2 = DelnucrCl2 + PartCl(I)*UncrossCl(I) DelnucrGa2 = DelnucrGa2 + PartGa(I)*UncrossGa(I) DellowerGa2 = DellowerGa2 + PartGa(I)*UnlowerGa(I) DelnucrLi2 = DelnucrLi2 + PartLi(I)*UncrossLi(I) End if End Do c c Add incoherently group 1 and 2. c DelnucrCl = (DelnucrCl1**2 + DelnucrCl2**2) DelnucrCl = sqrt(DelnucrCl) DelnucrGa = (DelnucrGa1**2 + DelnucrGa2**2) DelnucrGa = sqrt(DelnucrGa) DellowerGa = (DellowerGa1**2 + DellowerGa2**2) DellowerGa = sqrt(DellowerGa) DelnucrLi = (DelnucrLi1**2 + DelnucrLi2**2) DelnucrLi = sqrt(DelnucrLi) c c************************************* c Uncertainty in Be7-e-capture rate c************************************* c c Gruzinov and Bahcall ApJ, 490, 437 (1997) give a 1-sigma uncertainty c for the rate. c Be7eCl = UncertBe7e*PartCl(4) Be7eGa = UncertBe7e*PartGa(4) Be7eLi = UncertBe7e*PartLi(4) c c End of the dissolution and reassembly of the nuclear cross c section uncertainties. c c ****************************************************************** c Add quadratically the uncertainties. c ****************************************************************** c c**************************** c First the Cl uncertainties. c**************************** c x11Cl = Del11Cl*Del11Cl x33Cl = Del33Cl*Del33Cl x34Cl = Del34Cl*Del34Cl x114Cl = Del114Cl*Del114Cl x17Cl = Del17Cl*Del17Cl xhepCl = DelhepCl*DelhepCl xZXCl = DelZXCl*DelZXCl xlumCl = DellumCl*DellumCl xAgeCl = DelAgeCl*DelAgeCl xopcCl = DelopcCl*DelopcCl xnucrossCl = DelnucrCl*DelnucrCl xBe7eCl = Be7eCl*Be7eCl xdiffuseCl = DeldiffuseCl*DeldiffuseCl c c Perform quadratic addition of uncertainties. c UnCl2=x11Cl+x33Cl+x34Cl+x17Cl+xZXCl+xnucrossCl+xopcCl+x114Cl UnCl2 = UnCl2+ xBe7eCl + xlumCl + xageCl +xdiffuseCl + xhepCl c c************************************* c Now calculate Gallium uncertainties. c************************************* c x11Ga = Del11Ga*Del11Ga x33Ga = Del33Ga*Del33Ga x34Ga = Del34Ga*Del34Ga x114Ga = Del114Ga*Del114Ga x17Ga = Del17Ga*Del17Ga xhepGa = DelhepGa*DelhepGa xZXGa = DelZXGa*DelZXGa xlumGa = DellumGa*DellumGa xAgeGa = DelAgeGa*DelAgeGa xopcGa = DelopcGa*DelopcGa xnucrossGa = DelnucrGa*DelnucrGa xBe7eGa = Be7eGa*Be7eGa xdiffuseGa= DeldiffuseGa*DeldiffuseGa c c Perform quadratic addition of Ga uncertainties. c UnGa2=x11Ga+x33Ga+x34Ga+x17Ga+xZXGa+xnucrossGa+xopcGa+x114Ga UnGa2 = UnGa2+ xBe7eGa + xlumGa + xageGa +xdiffuseGa +xhepGa UnGalower2 = UnGa2 -xnucrossGa + DellowerGa*DellowerGa c c UnGalower2 uses the 3 sigma lower limit cross sections for gallium instead c of the somewhat larger upper limits c c******************************** c Now calculate Li uncertainties. c******************************** c x11Li = Del11Li*Del11Li x33Li = Del33Li*Del33Li x34Li = Del34Li*Del34Li x114Li = Del114Li*Del114Li x17Li = Del17Li*Del17Li xhepLi = DelhepLi*DelhepLi xZXLi = DelZXLi*DelZXLi xlumLi = DellumLi*DellumLi xAgeLi = DelAgeLi*DelAgeLi xopcLi = DelopcLi*DelopcLi xnucrossLi = DelnucrLi*DelnucrLi xBe7eLi = Be7eLi*Be7eLi xdiffuseLi= DeldiffuseLi*DeldiffuseLi c c Perform quadratic addition of Li uncertainties. c UnLi2=x11Li+x33Li+x34Li+x17Li+xZXLi+xnucrossLi+xopcLi+x114Li UnLi2 = UnLi2+ xBe7eLi + xlumLi + xageLi +xdiffuseLi +xhepLi c c****************************************************** c Calculate square root of the quadratic uncertainties. c***************************************************** c UnCl = sqrt(unCl2) UnGa = sqrt(unGa2) UnGalower = sqrt(UnGalower2) UnLi = sqrt(unLi2) c c ************************************** c Convert to 1 sigma instead of 3 sigma and insure that numbers are positive. c Then sum for the calculation of the average uncertainty. c ************************************** c c First calculate for chlorine c Del11Cl = abs(Del11Cl)/3.0 Del33Cl = abs(Del33Cl)/3.0 Del34Cl = abs(Del34Cl)/3.0 Del114Cl = abs(Del114Cl)/3.0 Del17Cl = abs(Del17Cl)/3.0 DelhepCl = abs(DelhepCl)/3.0 DelZXCl = abs(DelZXCl)/3.0 DellumCl = abs(DellumCl)/3.0 DelageCl = abs(DelageCl)/3.0 DelopcCl = abs(DelopcCl)/3.0 DeldiffuseCl = abs(DeldiffuseCl)/3.0 Be7eCl = abs(Be7eCl)/3.0 DelnucrCl = abs(DelnucrCl)/3.0 OneSigmaCl = UnCl/3.0 c c Calculate average uncertainties across all cases c K = CL av11(K) = av11(K) + Del11Cl av33(K) = av33(K) + Del33Cl av34(K) = av34(K) + Del34Cl av114(K) = av114(K) + Del114Cl av17(K) = av17(K) + Del17Cl avhep(K) = avhep(K) + DelhepCl avZX(K) = avZX(K) + DelZXCl avopac(K) = avopac(K) + DelopcCl avlum(K) = avlum(K) + DellumCl avage(K) = avage(K) + DelageCl avdiffuse(K) = avdiffuse(K) + DeldiffuseCl avUncert(K) = avUncert(K) + OneSigmaCl avBe7e(K) = avBe7e(K) + Be7eCl avnucr(K) = avnucr(K) + DelnucrCl c c Now calculate for gallium c Del11Ga = abs(Del11Ga)/3.0 Del33Ga = abs(Del33Ga)/3.0 Del34Ga = abs(Del34Ga)/3.0 Del114Ga = abs(Del114Ga)/3.0 Del17Ga = abs(Del17Ga)/3.0 DelhepGa = abs(DelhepGa)/3.0 DelZXGa = abs(DelZXGa)/3.0 DellumGa = abs(DellumGa)/3.0 DelageGa = abs(DelageGa)/3.0 DelopcGa = abs(DelopcGa)/3.0 DeldiffuseGa = abs(DeldiffuseGa)/3.0 Be7eGa = abs(Be7eGa)/3.0 DelnucrGa = abs(DelnucrGa)/3.0 OneSigmaGa = UnGa/3.0 UnGalower = UnGalower/3.0 c c Calculate average uncertainties across all cases c K = GA av11(K) = av11(K) + Del11Ga av33(K) = av33(K) + Del33Ga av34(K) = av34(K) + Del34Ga av114(K) = av114(K) + Del114Ga av17(K) = av17(K) + Del17Ga avhep(K) = avhep(K) + DelhepGa avZX(K) = avZX(K) + DelZXGa avopac(K) = avopac(K) + DelopcGa avlum(K) = avlum(K) + DellumGa avage(K) = avage(K) + DelageGa avdiffuse(K) = avdiffuse(K) + DeldiffuseGa avUncert(K) = avUncert(K) + OneSigmaGa avBe7e(K) = avBe7e(K) + Be7eGa avnucr(K) = avnucr(K) + DelnucrGa c c Convert Li to 1 sigma c Del11Li = abs(Del11Li)/3.0 Del33Li = abs(Del33Li)/3.0 Del34Li = abs(Del34Li)/3.0 Del114Li = abs(Del114Li)/3.0 Del17Li = abs(Del17Li)/3.0 DelhepLi = abs(DelhepLi)/3.0 DelZXLi = abs(DelZXLi) DelZXLi = abs(DelZXLi)/3.0 DellumLi = abs(DellumLi)/3.0 DelageLi = abs(DelageLi)/3.0 DelopcLi = abs(DelopcLi)/3.0 DeldiffuseLi = abs(DeldiffuseLi)/3.0 DelnucrLi = abs(DelnucrLi)/3.0 Be7eLi = abs(Be7eLi)/3.0 OneSigmaLi = UnLi/3.0 c c Calculate average uncertainties across all cases c K = LI av11(K) = av11(K) + Del11Li av33(K) = av33(K) + Del33Li av34(K) = av34(K) + Del34Li av114(K) = av114(K) + Del114Li av17(K) = av17(K) + Del17Li avhep(K) = avhep(K) + DelhepLi avZX(K) = avZX(K) + DelZXLi avopac(K) = avopac(K) + DelopcLi avlum(K) = avlum(K) + DellumLi avage(K) = avage(K) + DelageLi avdiffuse(K) = avdiffuse(K) + DeldiffuseLi avUncert(K) = avUncert(K) + OneSigmaLi avBe7e(K) = avBe7e(K) + Be7eLi avnucr(K) = avnucr(K) + DelnucrLi c c ******************************** c Uncertainties in neutrino fluxes c ******************************** c c Use improved version equation (24) of RMP 82. The version used here c reduces to the previous version when the individual uncertainties are c small. It is equation (16) of the Bahcall-Ulrich 1988 RMP paper. c c For all except B8, the uncertainty in the 7Be flux from the 7Be electron c capture rate is negligible because the 7Be electron capture is so fast, c about 10^3 times faster than the proton capture rate. c c Compute the sum of the squares and then take square root and divide c by three to get the 1 sigma uncertainties. c c Use a parameter to control the index for partial derivatives c c pp neutrino flux c K = PP c c Calculate individual contributions to pp uncertainty. c pp11 = (Uncert11**S11(K) - 1.0)*(Uncert11**S11(K) - 1.0) pp33 = (Uncert33**S33(K) - 1.0)*(Uncert33**S33(K) - 1.0) pp34 = (Uncert34**S34(K) - 1.0)*(Uncert34**S34(K) - 1.0) pp114 = (Uncert114**S114(K) - 1.0)*(Uncert114**S114(K) - 1.0) pp17 = (Uncert17**S17(K) - 1.0)*(Uncert17**S17(K) - 1.0) If (IZ .eq. 7) Then ppZX = UncertComposition(K)*UncertComposition(K) Else ppZX = (UncertZX**ZoverX(K) - 1.0)*(UncertZX**ZoverX(K) - 1.0) Endif ppopac = opc(K)*opc(K) pplum=(Uncertlum**sunlum(K)-1.0)*(Uncertlum**sunlum(K)-1.0) ppage = (Uncertage**age(K)-1.0)*(Uncertage**age(K)-1.0) ppdiffuse = Uncertdiffuse*Uncertdiffuse*diffuse(K)*diffuse(K) c c Sum the squares c Uncertsq(K) = pp11 + pp33 + pp34 + pp114 + pp17 +ppZX + * ppopac + pplum + ppage + ppdiffuse c c Take the square roots and adjust to 1 sigma. c Uncert(K) = sqrt(uncertsq(K))/3.0 pp11 = Sqrt(pp11)/3.0 pp33 = Sqrt(pp33)/3.0 pp34 = Sqrt(pp34)/3.0 pp114 = Sqrt(pp114)/3.0 pp17 = Sqrt(pp17)/3.0 ppZX = Sqrt(ppZX)/3.0 ppopac = Sqrt(ppopac)/3.0 pplum = sqrt(pplum)/3.0 ppage = sqrt(ppage)/3.0 ppdiffuse = sqrt(ppdiffuse)/3.0 c c Calculate average uncertainties across all cases c av11(K) = av11(K) + pp11 av33(K) = av33(K) + pp33 av34(K) = av34(K) + pp34 av114(K) = av114(K) + pp114 av17(K) = av17(K) + pp17 avZX(K) = avZX(K) + ppZX avopac(K) = avopac(K) + ppopac avlum(K) = avlum(K) + pplum avage(K) = avage(K) + ppage avdiffuse(K) = avdiffuse(K) + ppdiffuse avUncert(K) = avUncert(K) + Uncert(K) c c pep neutrino flux c K = PEP c c Calculate individual contributions to pep uncertainty. c pep11 = (Uncert11**S11(K) - 1.0)*(Uncert11**S11(K) - 1.0) pep33 = (Uncert33**S33(K) - 1.0)*(Uncert33**S33(K) - 1.0) pep34 = (Uncert34**S34(K) - 1.0)*(Uncert34**S34(K) - 1.0) pep114 = (Uncert114**S114(K) - 1.0)*(Uncert114**S114(K) - 1.0) pep17 = (Uncert17**S17(K) - 1.0)*(Uncert17**S17(K) - 1.0) If (IZ .eq. 7) Then pepZX = UncertComposition(K)*UncertComposition(K) Else pepZX = (UncertZX**ZoverX(K) - 1.0)*(UncertZX**ZoverX(K) - 1.0) Endif pepopac = opc(K)*opc(K) peplum=(Uncertlum**sunlum(K)-1.0)*(Uncertlum**sunlum(K)-1.0) pepage = (Uncertage**age(K)-1.0)*(Uncertage**age(K)-1.0) pepdiffuse = Uncertdiffuse*Uncertdiffuse*diffuse(K)*diffuse(K) c c Sum the squares c Uncertsq(K) = pep11 + pep33 + pep34 + pep114 + pep17 +pepZX + * pepopac + peplum + pepage + pepdiffuse c c Take the square roots and adjust to 1 sigma. c Uncert(K) = sqrt(uncertsq(K))/3.0 pep11 = Sqrt(pep11)/3.0 pep33 = Sqrt(pep33)/3.0 pep34 = Sqrt(pep34)/3.0 pep114 = Sqrt(pep114)/3.0 pep17 = Sqrt(pep17)/3.0 pepZX = Sqrt(pepZX)/3.0 pepopac = Sqrt(pepopac)/3.0 peplum = sqrt(peplum)/3.0 pepage = sqrt(pepage)/3.0 pepdiffuse = sqrt(pepdiffuse)/3.0 c c Calculate average uncertainties across all cases c av11(K) = av11(K) + pep11 av33(K) = av33(K) + pep33 av34(K) = av34(K) + pep34 av114(K) = av114(K) + pep114 av17(K) = av17(K) + pep17 avZX(K) = avZX(K) + pepZX avopac(K) = avopac(K) + pepopac avlum(K) = avlum(K) + peplum avage(K) = avage(K) + pepage avdiffuse(K) = avdiffuse(K) + pepdiffuse avUncert(K) = avUncert(K) + Uncert(K) c c Be7 neutrino flux c K = BE7 c c Calculate individual contributions to Be7 uncertainty. c Be711 = (Uncert11**S11(K) - 1.0)*(Uncert11**S11(K) - 1.0) Be733 = (Uncert33**S33(K) - 1.0)*(Uncert33**S33(K) - 1.0) Be734 = (Uncert34**S34(K) - 1.0)*(Uncert34**S34(K) - 1.0) Be7114 = (Uncert114**S114(K) - 1.0)*(Uncert114**S114(K) - 1.0) Be717 = (Uncert17**S17(K) - 1.0)*(Uncert17**S17(K) - 1.0) If (IZ .eq. 7) Then Be7ZX = UncertComposition(K)*UncertComposition(K) Else Be7ZX = (UncertZX**ZoverX(K) - 1.0)*(UncertZX**ZoverX(K) - 1.0) Endif Be7opac = opc(K)*opc(K) Be7lum=(Uncertlum**sunlum(K)-1.0)*(Uncertlum**sunlum(K)-1.0) Be7age = (Uncertage**age(K)-1.0)*(Uncertage**age(K)-1.0) Be7diffuse = Uncertdiffuse*Uncertdiffuse*diffuse(K)*diffuse(K) c c Sum the squares c Uncertsq(K) = Be711 + Be733 + Be734 + Be7114 + Be717 +Be7ZX + * Be7opac + Be7lum + Be7age + Be7diffuse c c Take the square roots and adjust to 1 sigma. c Uncert(K) = sqrt(uncertsq(K))/3.0 Be711 = Sqrt(Be711)/3.0 Be733 = Sqrt(Be733)/3.0 Be734 = Sqrt(Be734)/3.0 Be7114 = Sqrt(Be7114)/3.0 Be717 = Sqrt(Be717)/3.0 Be7ZX = Sqrt(Be7ZX)/3.0 Be7opac = Sqrt(Be7opac)/3.0 Be7lum = sqrt(Be7lum)/3.0 Be7age = sqrt(Be7age)/3.0 Be7diffuse = sqrt(Be7diffuse)/3.0 c c Calculate average uncertainties across all cases c av11(K) = av11(K) + Be711 av33(K) = av33(K) + Be733 av34(K) = av34(K) + Be734 av114(K) = av114(K) + Be7114 av17(K) = av17(K) + Be717 avZX(K) = avZX(K) + Be7ZX avopac(K) = avopac(K) + Be7opac avlum(K) = avlum(K) + Be7lum avage(K) = avage(K) + Be7age avdiffuse(K) = avdiffuse(K) + Be7diffuse avUncert(K) = avUncert(K) + Uncert(K) c c B8 neutrino flux c K = B8 c c Calculate individual contributions to B8 uncertainty. c B811 = (Uncert11**S11(K) - 1.0)*(Uncert11**S11(K) - 1.0) B833 = (Uncert33**S33(K) - 1.0)*(Uncert33**S33(K) - 1.0) B834 = (Uncert34**S34(K) - 1.0)*(Uncert34**S34(K) - 1.0) B8114 = (Uncert114**S114(K) - 1.0)*(Uncert114**S114(K) - 1.0) B817 = (Uncert17**S17(K) - 1.0)*(Uncert17**S17(K) - 1.0) If (IZ .eq. 7) Then B8ZX = UncertComposition(K)*UncertComposition(K) Else B8ZX = (UncertZX**ZoverX(K) - 1.0)*(UncertZX**ZoverX(K) - 1.0) Endif B8opac = opc(K)*opc(K) B8lum=(Uncertlum**sunlum(K)-1.0)*(Uncertlum**sunlum(K)-1.0) B8age = (Uncertage**age(K)-1.0)*(Uncertage**age(K)-1.0) B8diffuse = Uncertdiffuse*Uncertdiffuse*diffuse(K)*diffuse(K) c c Sum the squares c Uncertsq(K) = B811 + B833 + B834 + B8114 + B817 +B8ZX + * B8opac + B8lum + B8age + B8diffuse + * UncertBe7e * UncertBe7e c c UncertBe7e only affects the 8B neutrino flux. See notes at the c beginning of the code. c c Take the square roots and adjust to 1 sigma. c Uncert(K) = sqrt(uncertsq(K))/3.0 B811 = Sqrt(B811)/3.0 B833 = Sqrt(B833)/3.0 B834 = Sqrt(B834)/3.0 B8114 = Sqrt(B8114)/3.0 B817 = Sqrt(B817)/3.0 B8ZX = Sqrt(B8ZX)/3.0 B8opac = Sqrt(B8opac)/3.0 B8lum = sqrt(B8lum)/3.0 B8age = sqrt(B8age)/3.0 B8diffuse = sqrt(B8diffuse)/3.0 c c Calculate average uncertainties across all cases c av11(K) = av11(K) + B811 av33(K) = av33(K) + B833 av34(K) = av34(K) + B834 av114(K) = av114(K) + B8114 av17(K) = av17(K) + B817 avZX(K) = avZX(K) + B8ZX avopac(K) = avopac(K) + B8opac avlum(K) = avlum(K) + B8lum avage(K) = avage(K) + B8age avdiffuse(K) = avdiffuse(K) + B8diffuse avUncert(K) = avUncert(K) + Uncert(K) c c Compute average for UncertBe7e here, just for B8. Adjust to 1 sigma first. c avBe7e(K) = avBe7e(K) + Uncertbe7e/3.0 c c N13 neutrino flux c K = N13 c c Calculate individual contributions to N13 uncertainty. c N1311 = (Uncert11**S11(K) - 1.0)*(Uncert11**S11(K) - 1.0) N1333 = (Uncert33**S33(K) - 1.0)*(Uncert33**S33(K) - 1.0) N1334 = (Uncert34**S34(K) - 1.0)*(Uncert34**S34(K) - 1.0) N13114 = (Uncert114**S114(K) - 1.0)*(Uncert114**S114(K) - 1.0) N1317 = (Uncert17**S17(K) - 1.0)*(Uncert17**S17(K) - 1.0) If (IZ .eq. 7) Then N13ZX = UncertComposition(K)*UncertComposition(K) Else N13ZX = (UncertZX**ZoverX(K) - 1.0)*(UncertZX**ZoverX(K) - 1.0) Endif N13opac = opc(K)*opc(K) N13lum=(Uncertlum**sunlum(K)-1.0)*(Uncertlum**sunlum(K)-1.0) N13age = (Uncertage**age(K)-1.0)*(Uncertage**age(K)-1.0) N13diffuse = Uncertdiffuse*Uncertdiffuse*diffuse(K)*diffuse(K) c c Sum the squares c Uncertsq(K) = N1311 + N1333 + N1334 + N13114 + N1317 +N13ZX + * N13opac + N13lum + N13age + N13diffuse c c Take the square roots and adjust to 1 sigma. c Uncert(K) = sqrt(uncertsq(K))/3.0 N1311 = Sqrt(N1311)/3.0 N1333 = Sqrt(N1333)/3.0 N1334 = Sqrt(N1334)/3.0 N13114 = Sqrt(N13114)/3.0 N1317 = Sqrt(N1317)/3.0 N13ZX = Sqrt(N13ZX)/3.0 N13opac = Sqrt(N13opac)/3.0 N13lum = sqrt(N13lum)/3.0 N13age = sqrt(N13age)/3.0 N13diffuse = sqrt(N13diffuse)/3.0 c c Calculate average uncertainties across all cases c av11(K) = av11(K) + N1311 av33(K) = av33(K) + N1333 av34(K) = av34(K) + N1334 av114(K) = av114(K) + N13114 av17(K) = av17(K) + N1317 avZX(K) = avZX(K) + N13ZX avopac(K) = avopac(K) + N13opac avlum(K) = avlum(K) + N13lum avage(K) = avage(K) + N13age avdiffuse(K) = avdiffuse(K) + N13diffuse avUncert(K) = avUncert(K) + Uncert(K) c c O15 neutrino flux c K = O15 c c Calculate individual contributions to O15 uncertainty. c O1511 = (Uncert11**S11(K) - 1.0)*(Uncert11**S11(K) - 1.0) O1533 = (Uncert33**S33(K) - 1.0)*(Uncert33**S33(K) - 1.0) O1534 = (Uncert34**S34(K) - 1.0)*(Uncert34**S34(K) - 1.0) O15114 = (Uncert114**S114(K) - 1.0)*(Uncert114**S114(K) - 1.0) O1517 = (Uncert17**S17(K) - 1.0)*(Uncert17**S17(K) - 1.0) If (IZ .eq. 7) Then O15ZX = UncertComposition(K)*UncertComposition(K) Else O15ZX = (UncertZX**ZoverX(K) - 1.0)*(UncertZX**ZoverX(K) - 1.0) Endif O15opac = opc(K)*opc(K) O15lum=(Uncertlum**sunlum(K)-1.0)*(Uncertlum**sunlum(K)-1.0) O15age = (Uncertage**age(K)-1.0)*(Uncertage**age(K)-1.0) O15diffuse = Uncertdiffuse*Uncertdiffuse*diffuse(K)*diffuse(K) c c Sum the squares c Uncertsq(K) = O1511 + O1533 + O1534 + O15114 + O1517 +O15ZX + * O15opac + O15lum + O15age + O15diffuse c c Take the square roots and adjust to 1 sigma. c Uncert(K) = sqrt(uncertsq(K))/3.0 O1511 = Sqrt(O1511)/3.0 O1533 = Sqrt(O1533)/3.0 O1534 = Sqrt(O1534)/3.0 O15114 = Sqrt(O15114)/3.0 O1517 = Sqrt(O1517)/3.0 O15ZX = Sqrt(O15ZX)/3.0 O15opac = Sqrt(O15opac)/3.0 O15lum = sqrt(O15lum)/3.0 O15age = sqrt(O15age)/3.0 O15diffuse = sqrt(O15diffuse)/3.0 c c Calculate average uncertainties across all cases c av11(K) = av11(K) + O1511 av33(K) = av33(K) + O1533 av34(K) = av34(K) + O1534 av114(K) = av114(K) + O15114 av17(K) = av17(K) + O1517 avZX(K) = avZX(K) + O15ZX avopac(K) = avopac(K) + O15opac avlum(K) = avlum(K) + O15lum avage(K) = avage(K) + O15age avdiffuse(K) = avdiffuse(K) + O15diffuse avUncert(K) = avUncert(K) + Uncert(K) c c F17 neutrino flux c K = F17 c c Calculate individual contributions to F17 uncertainty. c F1711 = (Uncert11**S11(K) - 1.0)*(Uncert11**S11(K) - 1.0) F1733 = (Uncert33**S33(K) - 1.0)*(Uncert33**S33(K) - 1.0) F1734 = (Uncert34**S34(K) - 1.0)*(Uncert34**S34(K) - 1.0) F17114 = (Uncert114**S114(K) - 1.0)*(Uncert114**S114(K) - 1.0) F1717 = (Uncert17**S17(K) - 1.0)*(Uncert17**S17(K) - 1.0) If (IZ .eq. 7) Then F17ZX = UncertComposition(K)*UncertComposition(K) Else F17ZX = (UncertZX**ZoverX(K) - 1.0)*(UncertZX**ZoverX(K) - 1.0) Endif F17opac = opc(K)*opc(K) F17lum=(Uncertlum**sunlum(K)-1.0)*(Uncertlum**sunlum(K)-1.0) F17age = (Uncertage**age(K)-1.0)*(Uncertage**age(K)-1.0) F17diffuse = Uncertdiffuse*Uncertdiffuse*diffuse(K)*diffuse(K) F17116 = 3*(1.7/9.4)*3*(1.7/9.4) c c F17116 is uncertainty due to 16O(p,\gamma)17F reaction. Used uncertainty c given by Adelberger et al. (1998). c c Sum the squares c Uncertsq(K) = F1711 + F1733 + F1734 + F17114 + F1717 +F17ZX + * F17opac + F17lum + F17age + F17diffuse + F17116 c c Take the square roots and adjust to 1 sigma. c Uncert(K) = sqrt(uncertsq(K))/3.0 F1711 = Sqrt(F1711)/3.0 F1733 = Sqrt(F1733)/3.0 F1734 = Sqrt(F1734)/3.0 F17114 = Sqrt(F17114)/3.0 F1717 = Sqrt(F1717)/3.0 F17ZX = Sqrt(F17ZX)/3.0 F17opac = Sqrt(F17opac)/3.0 F17lum = sqrt(F17lum)/3.0 F17age = sqrt(F17age)/3.0 F17diffuse = sqrt(F17diffuse)/3.0 c c Calculate average uncertainties across all cases c av11(K) = av11(K) + F1711 av33(K) = av33(K) + F1733 av34(K) = av34(K) + F1734 av114(K) = av114(K) + F17114 av17(K) = av17(K) + F1717 avZX(K) = avZX(K) + F17ZX avopac(K) = avopac(K) + F17opac avlum(K) = avlum(K) + F17lum avage(K) = avage(K) + F17age avdiffuse(K) = avdiffuse(K) + F17diffuse avUncert(K) = avUncert(K) + Uncert(K) c c hep neutrino flux c K = HEP c c Calculate individual contributions to hep uncertainty. c hep11 = (Uncert11**S11(K) - 1.0)*(Uncert11**S11(K) - 1.0) hep33 = (Uncert33**S33(K) - 1.0)*(Uncert33**S33(K) - 1.0) hep34 = (Uncert34**S34(K) - 1.0)*(Uncert34**S34(K) - 1.0) hep114 = (Uncert114**S114(K) - 1.0)*(Uncert114**S114(K) - 1.0) hep17 = (Uncert17**S17(K) - 1.0)*(Uncert17**S17(K) - 1.0) hephep = (Uncerthep**Shep(K) - 1.0)*(Uncerthep**Shep(K) - 1.0) If (IZ .eq. 7) Then hepZX = UncertComposition(K)*UncertComposition(K) Else hepZX = (UncertZX**ZoverX(K) - 1.0)*(UncertZX**ZoverX(K) - 1.0) Endif hepopac = opc(K)*opc(K) heplum=(Uncertlum**sunlum(K)-1.0)*(Uncertlum**sunlum(K)-1.0) hepage = (Uncertage**age(K)-1.0)*(Uncertage**age(K)-1.0) hepdiffuse = Uncertdiffuse*Uncertdiffuse*diffuse(K)*diffuse(K) c c Sum the squares c Uncertsq(K) = hep11 + hep33 + hep34 + hep114 + hep17 +hepZX + * hepopac + heplum + hepage + hepdiffuse + hephep c c Take the square roots and adjust to 1 sigma. c Uncert(K) = sqrt(uncertsq(K))/3.0 hep11 = Sqrt(hep11)/3.0 hep33 = Sqrt(hep33)/3.0 hep34 = Sqrt(hep34)/3.0 hep114 = Sqrt(hep114)/3.0 hep17 = Sqrt(hep17)/3.0 hephep = Sqrt(hephep)/3.0 hepZX = Sqrt(hepZX)/3.0 hepopac = Sqrt(hepopac)/3.0 heplum = sqrt(heplum)/3.0 hepage = sqrt(hepage)/3.0 hepdiffuse = sqrt(hepdiffuse)/3.0 c c Calculate average uncertainties across all cases c av11(K) = av11(K) + hep11 av33(K) = av33(K) + hep33 av34(K) = av34(K) + hep34 av114(K) = av114(K) + hep114 av17(K) = av17(K) + hep17 avhep(K) = avhep(K) + hephep avZX(K) = avZX(K) + hepZX avopac(K) = avopac(K) + hepopac avlum(K) = avlum(K) + heplum avage(K) = avage(K) + hepage avdiffuse(K) = avdiffuse(K) + hepdiffuse avUncert(K) = avUncert(K) + Uncert(K) c c Adjust any remaining 3 sigma quantities to 1 sigma before printing c Uncert11 = (Uncert11 - 1.0)/3.0 Uncert33 = (Uncert33 - 1.0)/3.0 Uncert34 = (Uncert34 - 1.0)/3.0 Uncert17 = (Uncert17 - 1.0)/3.0 Uncerthep = (Uncerthep -1.0)/3.0 UncertZX = (UncertZX - 1.0)/3.0 Uncert114 = (Uncert114 - 1.0)/3.0 Uncertlum = (Uncertlum - 1.0)/3.0 Uncertage = (Uncertage - 1.0)/3.0 Uncertdiffuse = Uncertdiffuse/3.0 UncertBe7e = UncertBe7e/3.0 DellowerGa = DellowerGa/3.0 c c Write rates and uncertainties for solar neutrino experiments c Write(1,120) Write(1,121)RateCl,OneSigmaCl Write(1,122)RateGa,OneSigmaGa,UnGalower c Write(1,122)RateGa,OneSigmaGa,UnGalower,DellowerGa Write(1,123)RateLi,OneSigmaLi c c Write table of uncertainties c Write(1,130) Write(1,131) Uncert11,Uncert33,Uncert34,Uncert17,UncertZX, * Uncert114, Uncerthep Write(1,171) Write(1,133)Source(PP),pp11,pp33,pp34,pp17,ppZX,ppopac,pp114 Write(1,133)Source(PEP),pep11,pep33,pep34,pep17,pepZX,pepopac, * pep114 Write(1,133)Source(BE7),Be711,Be733,Be734,Be717,Be7ZX,Be7opac, * Be7114 Write(1,133)Source(B8),B811,B833,B834,B817,B8ZX,B8opac,B8114 Write(1,133)Source(N13),N1311,N1333,N1334,N1317,N13ZX,N13opac, * N13114 Write(1,133)Source(O15),O1511,O1533,O1534,O1517,O15ZX,O15opac, * O15114 Write(1,133)Source(F17),F1711,F1733,F1734,F1717,F17ZX,F17opac, * F17114 Write(1,170)Source(HEP),hep11,hep33,hep34,hep17,hepZX,hepopac, * hep114,hephep Write(1,133) Source(CL),Del11Cl,Del33Cl,Del34Cl,Del17Cl,DelZXCl, * DelopcCl,Del114Cl Write(1,133) Source(GA),Del11Ga,Del33Ga,Del34Ga,Del17Ga,DelZXGa, * DelopcGa,Del114Ga Write(1,133) Source(LI),Del11Li,Del33Li,Del34Li,Del17Li,DelZXLi, * DelopcLi,Del114Li Write(1,140) Write(1,141) Uncertlum,Uncertage,Uncertdiffuse,UncertBe7e Write(1,171) Write(1,142)Source(PP),pplum,ppage,ppdiffuse,Uncert(PP) Write(1,142)Source(PEP),peplum,pepage,pepdiffuse,Uncert(PEP) Write(1,142)Source(BE7),Be7lum,Be7age,Be7diffuse,Uncert(BE7) Write(1,144)Source(B8),B8lum,B8age,B8diffuse,UncertBe7e,Uncert(B8) Write(1,142)Source(N13),N13lum,N13age,N13diffuse,Uncert(N13) Write(1,142)Source(O15),O15lum,O15age,O15diffuse,Uncert(O15) Write(1,142)Source(F17),F17lum,F17age,F17diffuse,Uncert(F17) Write(1,142)Source(HEP),heplum,hepage,hepdiffuse,Uncert(HEP) Write(1,135) Source(CL),DellumCl,DelAgeCl,DeldiffuseCl, * Be7eCl,DelnucrCl,OneSigmaCl Write(1,135) Source(GA),DellumGa,DelAgeGa,DeldiffuseGa, * Be7eGa,DelnucrGa,OneSigmaGa Write(1,135) Source(LI),DellumLi,DelAgeLi,DeldiffuseLi, * Be7eLi,DelnucrLi,OneSigmaLi c c c****************************************************** c End of do loop over different cases of Uncertainties. c****************************************************** c End do c c********************************* c Print out average Uncertainties. c********************************* c 999 Continue c c Note: c Uncomment all write(1,400) statements to get top of page functionality c c Write(1,400) Write(1,101) Ncase Write(1,150) Write(1,171) Do K = 1, 11 av11(K) = av11(K)/Ncase av33(K) = av33(K)/Ncase av34(K) = av34(K)/Ncase av114(K) = av114(K)/Ncase av17(K) = av17(K)/Ncase avhep(K) = avhep(K)/Ncase avZX(K) = avZX(K)/Ncase avopac(K) = avopac(K)/Ncase avlum(K) = avlum(K)/Ncase avage(K) = avage(K)/Ncase avdiffuse(K) = avdiffuse(K)/Ncase avBe7e(K) = avBe7e(K)/Ncase avnucr(K) = avnucr(K)/Ncase avUncert(K) = avUncert(K)/Ncase Write(1,133)Source(K),av11(K),av33(K),av34(K),av17(K), * avZX(K), avopac(K),av114(K),avhep(K) Enddo c c Second half of table. B8 gets UncertBe7e printed. Cl, Ga and Li get c Be7e and nucrs values. c Write(1,160) Write(1,171) Do K = 1, 11 If(K.eq.B8) then Write(1,144)Source(K),avlum(K),avage(K),avdiffuse(K), * avBe7e(K),avUncert(K) Elseif(K.ge.CL) then Write(1,135)Source(K),avlum(K),avage(K),avdiffuse(K), * avBe7e(K),avnucr(K),avUncert(K) Else Write(1,142)Source(K),avlum(K),avage(K),avdiffuse(K), * avUncert(K) Endif Enddo c c ******************************** c Format Statements c ******************************** c 100 Format(/' Solar Neutrino Fluxes, Rates, and Uncertainties') 101 Format(/' **********************'/ * ' Average over ', I2, ' cases.'/ * ' **********************') 102 Format(/' Case number: ',I2) 103 Format( *' SNU Contributions'// *' Source Fluxes Cl Ga Li'/) 104 Format(1x,A8,E10.3, F10.2, F10.1, F10.1) 120 Format(//' The Capture Rates and Uncertainties for Solar Neutrino * Experiments: ',/) 121 Format(' Cl rate = ',F7.2,' +/-',F5.2) c122 Format(' Ga rate = ',F7.2,' + ',F5.2,' -',F5.2, c * ' (LowerGanucr =',F5.2,')') 122 Format(' Ga rate = ',F7.2,' + ',F5.2,' -',F5.2) 123 Format(' Li rate = ',F7.2,' +/-',F5.2) 130 Format(// *' Fractional 1-Sigma Uncertainties from different input data'// *' Source 11 33 34 17 Z/X Opac 114 hep'/) 131 Format(' Uncert ',4F6.3,2x,F6.3,6x,F6.3,1x,F6.3/) 133 Format(1x, A8, 4F6.3,2x,4F6.3,5F6.3) 135 Format(1x, A8, 3F6.3,2x,2F6.3,1x,F6.3) 140 Format(// *' Fractional 1-Sigma Uncertainties from different input data'// *' Source Lumin Age Diff Be7e NuCrs Total'/) 141 Format(' Uncert ',4F6.3/) 142 Format(1x, A8, 3F6.3,15x,F6.3) 144 Format(1x, A8, 3F6.3,2x,F6.3,7x,F6.3) 150 Format(//' Averages:', *' Fractional 1-Sigma Uncertainties from different input data'// *' Source 11 33 34 17 Z/X Opac 114 hep'/) 160 Format(//' Averages:', *' Fractional 1-Sigma Uncertainties from different input data'// *' Source Lumin Age Diff Be7e NuCrs Total'/) 170 Format(1x, A8, 4F6.3,2x,3F6.3,2x,F6.3,6F6.3) 171 Format(1x, 'Fractional uncertainties in nu fluxes; Delta SNU for *detectors',/) 180 Format(1x,//' Neutrino Flux and Rate Uncertainties Computed Using * Individual Abundances or Total Z/X.',//) 181 Format(1x,/,'1sigma individual composition uncertainties and 1also radiochemical rate uncertainties.') 182 Format(1x,/,'Neutrino fluxes',//) 184 Format(1x,//,8F10.4,/, 3F10.2,//) 185 Format(1x,/, 'If IZ = 7,then individual compositions; IZ = ',I4,/) 300 Format(/////////A1) 400 Format('1') Stop End