|
GSI
FRS
| |
For Windows/XP users > Run Karol ( P:\frstools\karol.bat
)
Windows users outside GSI need the executable and some additional
files.
Fortran version of the code is
also available.
Source Text (PL/I)
Module Karol_Main
%PROCESS MACRO LIMITS(EXTNAME(31)) LANGLVL(SAA2) MARGINS(1,100);
Karol_Main: Procedure Options(Main);
/* Main Procedure to call the Karol Model for calculating */
/* the total nuclear cross sections */
Declare /**/
(I_Z1,I_Z2) Bin Fixed(16), /* Projectile */
(I_A1,I_A2) Bin Fixed(16), /* Target */
R_EA1 Dec Float(16),
/* Projectile energy per nucleon in MeV */
R_SigTot Dec Float(16),
C_Dummy Char(8) Var;
Declare /**/
Karol Entry(Bin Fixed(16), Bin Fixed(16),
Bin Fixed(16), Bin Fixed(16),
Dec Float(16)) Returns(Dec FLoat(16));
Put Edit('Calculation of the total nuclear interaction') (Skip,A);
Put Edit('cross section') (Skip,A);
Put Edit('according to Karol (Phys. Rev. C11 (1975) 1203)') (Skip,A);
Put Edit('with modifications introduced by Brohm') (Skip,A);
Put Edit('(Nucl. Phys. A 569 (1994) 821)') (Skip,A);
Put Edit('Enter Z and A of projectile: ') (Skip,Skip,A);
Get List(I_Z1,I_A1);
Put Edit('Enter Z and A of target: ') (A);
Get List(I_Z2,I_A2);
Put Edit('Enter projectile energy per nucleon in MeV: ') (A);
Get List(R_EA1);
R_SigTot = Karol(I_Z1,I_A1,I_Z2,I_A2,R_EA1);
Put Edit('The cross section (in mb) is ') (Skip,A);
Put Edit(R_SigTot) (F(10));
Get Edit(C_Dummy) (A(8));
End Karol_Main;
Module Karol
%PROCESS MACRO LIMITS(EXTNAME(31)) LANGLVL(SAA2) MARGINS(1,100);
O_Karol: Package Exports(Karol);
Declare /* Global values */
(R_ZP,R_AP) Dec Float(16) Static, /* Projectile */
(R_ZT,R_AT) Dec Float(16) Static, /* Target */
R_EAP Dec Float(16) Static, /* Energy per nucleon in MeV */
R_LimP Dec Float(16) Static, /* Limit for integration, proj. */
R_LimT Dec Float(16) Static, /* Limit for integration, target */
R_Eps Dec Float(16) Static, /* Required precision */
(TNP,TNT,TPP,TPT) Dec Float(16) Static,
(SIGPP,SIGNP) Dec Float(16) Static,
RXA Dec Float(16) Static,
RYA Dec Float(16) Static,
RBA Dec Float(16) Static,
DRHOPTA Dec Float(16) Static,
DRHONTA Dec Float(16) Static,
DRHOPPA Dec Float(16) Static,
DRHONPA Dec Float(16) Static,
LDMK Dec Float(16) Static,
LDMA Dec Float(16) Static,
LDMC Dec Float(16) Static,
LDMQ Dec Float(16) Static,
LDMJ Dec Float(16) Static,
LDML Dec Float(16) Static,
LDMR Dec Float(16) Static,
LDMD Dec Float(16) Static,
Pi Dec Float(16) Static Init(3.141592654E0);
Declare /**/
(Float,Exp) Builtin;
/*********************************************************************/
Karol: Procedure(I_ZP_in,I_AP_in,I_ZT_in,I_AT_in,R_EAP_in)
Returns(Dec Float(16));
Declare /**/
(I_ZP_in,I_AP_in) Bin Fixed(15), /* Projectile */
(I_ZT_in,I_AT_in) Bin Fixed(15), /* Target */
R_EAP_in Dec Float(16), /* Energy per nucleon in MeV */
R_Cutoff Dec Float(16), /* Cutoff parameter */
R_Bintm Dec Float(16), /* Max impact parameter for integr. */
R_CrNuc Dec Float(16); /* Interaction cross section */
/* TOTAL NN CROSS SECTION FROM: LOCK + MEASDAY, INTERMEDIATE ENERGY */
/* NUCLEAR PHYSICS, METHUEN, LONDON 1970, P. 180 FF. */
%DECLARE N1TAB FIXED;
%N1TAB = 12;
%DECLARE N2TAB FIXED;
%N2TAB = 12;
Declare /**/
E1TAB(N1TAB) Dec Float(16) Init(
10., 20., 50., 100., 150., 250., 300., 400., 500., 700., 800.,
1000.),
S1TAB(N1TAB) Dec Float(16) Init(
300., 150., 52., 28., 24., 23., 23.5, 25.5, 30.5, 44.5, 47.,48.),
E2TAB(N2TAB) Dec Float(16) Init(
10., 20., 50., 100., 150., 250., 300., 400., 500., 700., 800.,
1000. ),
S2TAB(N2TAB) Dec Float(16) Init(
1000., 500., 167., 77., 50., 36., 35., 33., 33.5,
36., 38., 40.);
/* Fill global values */
IF I_ZP_in RINT1DL: NEGATIVE VALUES NOT ALLOWED!') (SKIP,A);
Return(0);
END;
DO I = 1 TO N;
IF( XTAB(I) <0 ) | (YTAB(I) < 0 ) THEN DO; Put Edit(' RINT1DL: NEGATIVE VALUES NOT ALLOWED!')
(SKIP,A);
Return(0);
END;
IF( XTAB(I) = 0 ) Then XTAB(I) = 1.E-50;
IF( YTAB(I) = 0 ) Then YTAB(I) = 1.E-50;
IF( X = XTAB(I) ) Then Do;
Y = LOG(YTAB(I));
Return(Exp(y));
END;
END;
IF( X XTAB(N) ) THEN DO;
ISAVE = N - 1;
GOTO L100;
END;
DO I = 1 To N;
IF( X 30.E0 Then TRA = 0;
Else TRA = Exp(-RES);
R_Value = 10.E0 * 2.E0 * Pi * (1.E0 - TRA) * RBA;
Return(R_Value);
End U_SiYint;
/*********************************************************************/
U_SiXint: Procedure(RY) Returns(Dec Float(16));
/* Integration perpendicular to impact parameter */
Declare /**/
RY Dec Float(16),
R_Value Dec Float(16);
RYA = RY;
R_Value = U_Simpson(-R_LimP, 0.E0, R_eps, U_SI) * 2.E0;
Return(R_Value);
End U_SiXint;
/*********************************************************************/
U_SI: Procedure(RX) Returns(Dec Float(16));
/* Functions for total cross section */
Declare /**/
RX Dec Float(16),
DPROB Dec Float(16),
DRHOPT Dec Float(16),
DRHONT Dec Float(16),
DRHOPP Dec Float(16),
DRHONP Dec Float(16);
RXA = RX;
/* Integrate densities over Z for given RXA, RYA and RBA */
DRHOPT = U_Simpson(-R_LimT, 0.E0, R_eps, U_RHOPT) * 2.E0;
DRHOPTA = DRHOPT;
DRHONT = U_Simpson(-R_LimT, 0.E0, R_eps, U_RHONT) * 2.E0;
DRHONTA = DRHONT;
DRHOPP = U_Simpson(-R_LimP, 0.E0, R_eps, U_RHOPP) * 2.E0;
DRHOPPA = DRHOPP;
DRHONP = U_Simpson(-R_LimP, 0.E0, R_eps, U_RHONP) * 2.E0;
DRHONPA = DRHONP;
/* Probability for collisions in this row */
DPROB = (DRHOPP*DRHOPT+0.94E0*DRHONP*DRHONT)*SIGPP
+(DRHOPP*DRHONT+DRHONP*DRHOPT)*SIGNP;
Return(DPROB);
End U_SI;
/*********************************************************************/
/* NUCLEAR DENSITY DISTRIBUTIONS:
F. Vives, K.-H. Schmidt, 15. 9. 2000
For A<=20 Gaussian distributions, RMS experimental values, if available from H. de Vries et al. Atomic Data and Nuclear Data tables, Vol. 36, No 3, May 1987 FOR A<="10" VOLUME-NORMALIZED GAUSSIAN (OR SURFACE-NORMALIZED IF RMS RADIUS KNOWN) LITERATURE: KAROL, PHYS. REV. C11 (1975) 1203. FOR A>10 DROPLET-DENSITIES ACCORDING TO
MYERS AND SWIATECKI, ANN. OF. PHYS. 84 (1974) 186,
DROPLET PARAMETERS SEE BERDICHEVSKY AND TONDEUR,
Z. PHYS. A322 (1985) 141, AND DIFFUSENESS ACCORDING
TO MYERS AND GROOTE, PHYS. LETT. 61B (1976) 125 */
U_RHOPT: Procedure(ZST) Returns(Dec Float(16));
Declare /**/
ZST Dec Float(16),
RHOPT Dec Float(16),
R13 Dec Float(16) Init(0.333333333333333E0),
R23 Dec Float(16) Init(0.666666666666667E0),
R43 Dec Float(16) Init(1.333333333333333E0),
RG Dec Float(16),
BPT Dec Float(16),
EPST Dec Float(16),
DELTAT Dec Float(16),
RT Dec Float(16),
DT Dec Float(16),
RPT Dec Float(16),
CPT Dec Float(16),
RHO0PT Dec Float(16),
AUT Dec Float(16),
RMS(4,20) Dec Float(16) INIT(0.84, 1.9, 0., 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2.1, 1.67,
2.43, 0., 2.45, 2.41, 2.54, 2.72, 0., 0., 0., 2.95, 2.95,
3.04, 0., 3.21, 0., 0., 0., 3.5, 1.7, 0., 2.39, 2.5, 2.42,
2.44, 2.54, 0., 0., 0., 0., 0., 0., 0., 3.07, 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., 2.77, 0., 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 3.50);
/* First index: N-Z+2
Second index: Z */
/* Some examples: */
/* Z N A N-Z N-Z+2
1 0 1 -1 1 proton
1 1 2 0 2 deuteron
1 2 3 1 3 triton */
/* RMS is an array with the value of RMS corresponding to a certain
Z value (column) and N value (row), as an example the rms value
for a Deuterium projectile is 2.1 and for a Tritium is 1.7 */
/* RMS CHARGE RADII FORM LANDOLT-BOERNSTEIN, VOL. 2, P.30FF */
RG = SQRT(RXA**2+(RBA+RYA)**2+ZST**2);
IF( R_AT <21.e0 ) THEN DO; /* DO i="-1" to 2 by 1 ; */ IF( RMS(FIXED(R_AT 2.E0 * R_ZT + 2. ,15),FIXED(R_ZT ,15))> 1.E-1 ) THEN DO;
AUT = RMS(FIXED( R_AT - 2.E0 * R_ZT + 2. ,15),FIXED(R_ZT ,15)) * 8.165E-1;
RHO0PT = R_ZT / ( AUT * 1.7725E0)**3;
END;
ELSE DO;
CPT = 1.07E0 * R_AT**R13;
AUT = SQRT(( 4.0E0 * CPT * TPT + TPT**2 ) / 6.43775E0 );
RHO0PT = 0.5E0 * 3.E0 * R_ZT
/ ( 4.E0 * PI * CPT**3
* ( 1.E0 + ( PI**2 * TPT**2 / (19.36E0 * CPT**2) )))
* EXP( (CPT / AUT)**2);
END;
IF( RG/AUT >= 5 ) THEN RHOPT = 0.E0;
IF( RG/AUT <5 ) THEN RHOPT="RHO0PT" * EXP( (RG/AUT)**2 ); END; ELSE DO; BPT="0.413E0" * TPT; DELTAT="((R_AT-2.E0*R_ZT)/R_AT+" (3.E0/16.E0)*(LDMC/LDMQ)*R_ZT*R_AT**(-R23)) / ( 1.E0 + (9.E0/4.E0)*(LDMJ/LDMQ) * R_AT**(-R13) ); EPST="(" 2.E0*LDMA*R_AT**(-R13)+LDML*DELTAT**2 + LDMC*R_ZT**2*R_AT**(-R43) )/LDMK; RT="LDMR" * R_AT**R13 * ( 1.E0 + EPST ); DT="R23" * ((R_AT-2.E0*R_ZT)/R_AT DELTAT ) * RT; RPT="RT" ((R_AT-R_ZT)/R_AT) * DT; CPT="RPT" * ( 1.E0 (BPT/RPT)**2 ); RHO0PT="3.E0" * R_ZT / ( 4.E0*PI * CPT**3 * ( 1.E0 + PI**2 * TPT**2 / (19.36E0 * CPT**2))); IF( (RG CPT)> 80 ) THEN RHOPT = 0.E0;
ELSE DO;
RHOPT = RHO0PT / ( 1.E0 + EXP((RG-CPT)/(TPT/4.4E0)) );
/* RHOPT = RHO0PT / ( 1.E0 + (EXP((RG-CPT)/(TPT/4.4E0)))**1.5 );
*/
END;
END;
IF( RHOPT <1.e-30 ) THEN RHOPT="0.E0;" RETURN(RHOPT); End U_RHOPT; U_RHONT: Procedure(ZST) Returns(Dec Float(16)); Declare /**/ ZST Dec Float(16), RHONT Dec Float(16), R13 Dec Float(16) Init(0.333333333333333E0), R23 Dec Float(16) Init(0.666666666666667E0), R43 Dec Float(16) Init(1.333333333333333E0), RG Dec Float(16), BNT Dec Float(16), EPST Dec Float(16), DELTAT Dec Float(16), RT Dec Float(16), DT Dec Float(16), RNT Dec Float(16), CNT Dec Float(16), RHO0NT Dec Float(16), AUT Dec Float(16), RMS(4,20) Dec Float(16) INIT(0.84, 1.9, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2.1, 1.67, 2.43, 0., 2.45, 2.41, 2.54, 2.72, 0., 0., 0., 2.95, 2.95, 3.04, 0., 3.21, 0., 0., 0., 3.5, 1.7, 0., 2.39, 2.5, 2.42, 2.44, 2.54, 0., 0., 0., 0., 0., 0., 0., 3.07, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2.77, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 3.50); /* RMS is an array with the value of RMS corresponding to a certain Z value (column) and N value (row), as an example the rms value for a Deuterium projectile is 2.1 and for a Tritium is 1.7 */ /* RMS CHARGE RADII FORM LANDOLT-BOERNSTEIN, VOL. 2, P.30FF */ RG="SQRT(RXA**2+(RBA+RYA)**2+ZST**2);" IF( R_AT < 21.E0 ) THEN DO; /* DO i="-1" to 2; */ IF( RMS(FIXED(R_AT 2.E0 * R_ZT + 2.,15),FIXED(R_ZT ,15))> 1.E-1 ) THEN DO;
AUT = RMS(FIXED(R_AT - 2.E0 * R_ZT + 2. ,15),FIXED(R_ZT,15)) * 8.165E-1;
RHO0NT = (R_AT - R_ZT) / ( AUT * 1.7725E0)**3;
END;
ELSE DO;
CNT = 1.07E0 * R_AT**R13;
AUT = SQRT(( 4.0E0 * CNT * TNT + TNT**2 ) / 6.43775E0 );
RHO0NT = 0.5E0 * 3.E0 * (R_AT-R_ZT)
/ ( 4.E0 * PI * CNT**3
* ( 1.E0 + ( PI**2 * TNT**2 / (19.36E0 * CNT**2) )))
* EXP( (CNT / AUT)**2);
END;
IF( RG/AUT >= 5 ) THEN RHONT = 0.E0;
IF( RG/AUT <5 ) THEN RHONT="RHO0NT" * EXP( (RG/AUT)**2 ); END; ELSE DO; BNT="0.413E0" * TNT; DELTAT="((R_AT-2.E0*R_ZT)/R_AT+" (3.E0/16.E0)*(LDMC/LDMQ)*R_ZT*R_AT**(-R23)) / ( 1.E0 + (9.E0/4.E0)*(LDMJ/LDMQ) * R_AT**(-R13) ); EPST="(" 2.E0*LDMA*R_AT**(-R13)+LDML*DELTAT**2 + LDMC*R_ZT**2*R_AT**(-R43) )/LDMK; RT="LDMR" * R_AT**R13 * ( 1.E0 + EPST ); DT="R23" * ((R_AT-2.E0*R_ZT)/R_AT DELTAT ) * RT; RNT="RT" + (R_ZT/R_AT) * DT; CNT="RNT" * ( 1.E0 (BNT/RNT)**2 ); RHO0NT="3.E0" * ( R_AT R_ZT )/ ( 4.E0*PI * CNT**3 * ( 1.E0 + PI**2 * TNT**2 / (19.36E0 * CNT**2))); IF( (RG CNT)> 80. ) THEN RHONT = 0.E0;
ELSE DO;
RHONT = RHO0NT / ( 1.E0 + EXP((RG-CNT)/(TNT/4.4E0)) );
/* RHONT = RHO0NT / ( 1.E0 + (EXP((RG-CNT)/(TNT/4.4E0)))**1.5 );
*/
END;
END;
IF( RHONT <1.e-30 ) THEN RHONT="0.E0;" RETURN(RHONT); End U_RHONT; U_RHOPP: Procedure(ZST) Returns(Dec Float(16)); Declare /**/ ZST Dec Float(16), RHOPP Dec Float(16), R13 Dec Float(16) Init(0.333333333333333E0), R23 Dec Float(16) Init(0.666666666666667E0), R43 Dec Float(16) Init(1.333333333333333E0), RG Dec Float(16), BPP Dec Float(16), EPSP Dec Float(16), DELTAP Dec Float(16), RP Dec Float(16), DP Dec Float(16), RPP Dec Float(16), CPP Dec Float(16), RHO0PP Dec Float(16), AUP Dec Float(16), RMS(4,20) Dec Float(16) INIT(0.84, 1.9, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2.1, 1.67, 2.43, 0., 2.45, 2.41, 2.54, 2.72, 0., 0., 0., 2.95, 2.95, 3.04, 0., 3.21, 0., 0., 0., 3.5, 1.7, 0., 2.39, 2.5, 2.42, 2.44, 2.54, 0., 0., 0., 0., 0., 0., 0., 3.07, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2.77, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 3.50); /* RMS is an array with the value of RMS corresponding to a certain Z value (column) and N value (row), as an example the rms value for a Deuterium projectile is 2.1 and for a Tritium is 1.7 */ /* RMS CHARGE RADII FORM LANDOLT-BOERNSTEIN, VOL. 2, P.30FF */ RG="SQRT(RXA**2+RYA**2+ZST**2);" IF( R_AP < 21.E0 ) THEN DO; /* DO i="-1" to 2 by 1; */ IF( RMS(FIXED(R_AP 2.E0 * R_ZP + 2. ,15),FIXED(R_ZP ,15))> 1.E-1 ) THEN DO;
AUP = RMS(FIXED(R_AP - 2.E0 * R_ZP + 2. ,15),FIXED(R_ZP,15)) * 8.165E-1;
RHO0PP = R_ZP / ( AUP * 1.7725E0)**3;
END;
ELSE DO;
CPP = 1.07E0 * R_AP**R13;
AUP = SQRT(( 4.0E0 * CPP * TPP + TPP**2 ) / 6.43775E0 );
RHO0PP = 0.5E0 * 3.E0 * R_ZP
/ ( 4.E0 * PI * CPP**3
* ( 1.E0 + ( PI**2 * TPP**2 / (19.36E0 * CPP**2) )))
* EXP( (CPP / AUP)**2);
END;
IF( RG/AUP >= 5 ) THEN RHOPP = 0.E0;
IF( RG/AUP <5 ) THEN RHOPP="RHO0PP" * EXP( (RG/AUP)**2 ); END; ELSE DO; BPP="0.413E0" * TPP; DELTAP="((R_AP-2.E0*R_ZP)/R_AP+" (3.E0/16.E0)*(LDMC/LDMQ)*R_ZP*R_AP**(-R23)) / ( 1.E0 + (9.E0/4.E0)*(LDMJ/LDMQ) * R_AP**(-R13) ); EPSP="(" 2.E0*LDMA*R_AP**(-R13)+LDML*DELTAP**2 + LDMC*R_ZP**2*R_AP**(-R43) )/LDMK; RP="LDMR" * R_AP**R13 * ( 1.E0 + EPSP ); DP="R23" * ((R_AP-2.E0*R_ZP)/R_AP DELTAP ) * RP; RPP="RP" ((R_AP-R_ZP)/R_AP) * DP; CPP="RPP" * ( 1.E0 (BPP/RPP)**2 ); RHO0PP="3.E0" * R_ZT / ( 4.E0*PI * CPP**3 * ( 1.E0 + PI**2 * TPP**2 / (19.36E0 * CPP**2))); IF( (RG CPP)> 80. ) THEN RHOPP = 0.E0;
ELSE DO;
RHOPP = RHO0PP / ( 1.E0 + EXP((RG-CPP)/(TPP/4.4E0)) );
/* RHOPP = RHO0PP / ( 1.E0 + (EXP((RG-CPP)/(TPP/4.4E0)))**1.5 );
*/
END;
END;
IF( RHOPP <1.e-30 ) THEN RHOPP="0.E0;" RETURN(RHOPP); End U_RHOPP; U_RHONP: Procedure(ZST) Returns(Dec Float(16)); Declare /**/ ZST Dec Float(16), RHONP Dec Float(16), R13 Dec Float(16) Init(0.333333333333333E0), R23 Dec Float(16) Init(0.666666666666667E0), R43 Dec Float(16) Init(1.333333333333333E0), RG Dec Float(16), BNP Dec Float(16), EPSP Dec Float(16), DELTAP Dec Float(16), RP Dec Float(16), DP Dec Float(16), RNP Dec Float(16), CNP Dec Float(16), RHO0NP Dec Float(16), AUP Dec Float(16), RMS(4,20) Dec Float(16) INIT(0.84, 1.9, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2.1, 1.67, 2.43, 0., 2.45, 2.41, 2.54, 2.72, 0., 0., 0., 2.95, 2.95, 3.04, 0., 3.21, 0., 0., 0., 3.5, 1.7, 0., 2.39, 2.5, 2.42, 2.44, 2.54, 0., 0., 0., 0., 0., 0., 0., 3.07, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 2.77, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 3.50); /* RMS is an array with the value of RMS corresponding to a certain Z value (row) and N value (column), as an example the rms value for a Deuterium projectile is 2.1 and for a Tritium is 1.7 */ /* RMS CHARGE RADII FORM LANDOLT-BOERNSTEIN, VOL. 2, P.30FF */ RG="SQRT(RXA**2+RYA**2+ZST**2);" IF( R_AP < 21.E0 ) THEN DO; /* DO i="-1" to 2 by 1; */ IF( RMS(FIXED(R_AP 2.E0 * R_ZP + 2.,15),FIXED(R_ZP ,15))> 1.E-1 ) THEN DO;
AUP = RMS(FIXED(R_AP - 2.E0 * R_ZP + 2. ,15),FIXED(R_ZP ,15)) * 8.165E-1;
RHO0NP =(R_AP - R_ZP) / ( AUP * 1.7725E0)**3;
END;
ELSE DO;
CNP = 1.07E0 * R_AP**R13;
AUP = SQRT(( 4.0E0 * CNP * TNP + TNP**2 ) / 6.43775E0 );
RHO0NP = 0.5E0 * 3.E0 * (R_AP - R_ZP)
/ ( 4.E0 * PI * CNP**3
* ( 1.E0 + ( PI**2 * TNP**2 / (19.36E0 * CNP**2) )))
* EXP( (CNP / AUP)**2);
END;
IF( RG/AUP >= 5 ) THEN RHONP = 0.E0;
IF( RG/AUP <5 ) THEN RHONP="RHO0NP" * EXP( (RG/AUP)**2 ); END; ELSE DO; BNP="0.413E0" * TNP; DELTAP="((R_AP-2.E0*R_ZP)/R_AP+" (3.E0/16.E0)*(LDMC/LDMQ)*R_ZP*R_AP**(-R23)) / ( 1.E0 + (9.E0/4.E0)*(LDMJ/LDMQ) * R_AP**(-R13) ); EPSP="(" 2.E0*LDMA*R_AP**(-R13)+LDML*DELTAP**2 + LDMC*R_ZP**2*R_AP**(-R43) )/LDMK; RP="LDMR" * R_AP**R13 * ( 1.E0 + EPSP ); DP="R23" * ((R_AP-2.E0*R_ZP)/R_AP DELTAP ) * RP; RNP="RP" + (R_ZP/R_AP) * DP; CNP="RNP" * ( 1.E0 (BNP/RNP)**2 ); RHO0NP="3.E0" * (R_AP R_ZT) / ( 4.E0*PI * CNP**3 * ( 1.E0 + PI**2 * TNP**2 / (19.36E0 * CNP**2))); IF( (RG CNP)> 80. ) THEN RHONP = 0.E0;
ELSE DO;
RHONP = RHO0NP / ( 1.E0 + EXP((RG-CNP)/(TNP/4.4E0)) );
/* RHONP = RHO0NP / ( 1.E0 + (EXP((RG-CNP)/(TNP/4.4E0)))**1.5 );
*/
END;
END;
IF( RHONP <1.e-30 ) THEN RHONP="0.E0;" RETURN(RHONP); End U_RHONP; /*********************************************************************/ U_Simpson: Procedure(A,B,EPS,FUNCTION) Returns(Dec Float(16)) ; /* INTEGRATES NUMERICALLY POSITIVE FUNCTION */ /* 0 <="A" <="X" <="B" */ /* KLAUS PIELENZ 1977 */ /* DECLARATION: */ /* DCL SIMPSON ENTRY(Dec Float(16), Dec Float(16), */ /* Dec Float(16),ENTRY) RETURNS(Dec Float(16)); */ DCL /**/ (A,B,EPS) Dec Float(16), /* PARAMETERS : */ /* INTEGRATION LIMITS, */ /* RELATIVE ERROR OF RESULT AND */ /* NAME OF FUNCTION */ FUNCTION ENTRY(Dec Float(16)) RETURNS(Dec Float(16)) VARIABLE, (H,S1,S2,S4,S,S_LAST) Dec Float(16), /* INTERVAL-WIDTH, HELP VARIABLES FOR ITERATION, */ /* TEMP. & LAST VALUE OF INTEGRAL */ K_S Bin Fixed(15), L Bin Fixed(31), /* LOCAL LOOP COUNTERS */ ABS BUILTIN; H="0.5E0" * (B A); if H <="0" then return(0); S1="H/3.E0" * ( FUNCTION(A) + FUNCTION(B) ); S2="0;" S4="4.E0*H/3.E0" * FUNCTION(A+H); S="S1" + S2 + S4; S_LAST="0;" /* READY FOR ITERATION */ Iter: DO L="4" REPEAT 2 * L WHILE(ABS(S-S_LAST)> EPS*S);
/* ITERATION-LOOP */
H = 0.5E0 * H;
/* HALF INTERVAL-WIDTH */
S1 = S1 * 0.5E0;
S2 = S2 * 0.5E0 + S4 * 0.25E0;
S4 = 0;
/* RESET FOR K_S-LOOP */
DO K_S=1 BY 2 TO L-1;
S4 = S4 + 4.E0*H/3.E0 * FUNCTION(A + K_S*H);
END;
S_LAST = S;
/* UPDATING */
S = S1 + S2 + S4;
IF L >= 2**30 THEN DO;
/* AVOID FIXEDOVERFLOW(L) */
Put edit('SPECIFIED ACCURACY NOT AVAILABLE') (SKIP,A);
LEAVE Iter;
END;
END;
RETURN(S);
END U_Simpson;
/*********************************************************************/
End O_Karol;
|
