CC CC ! ! CC ! G E K Ü R Z T E U N D V E R Ä N D E R T E V E R S I O N ! CC ! ! CC ! F Ü R M C V 3 K ! CC ! ! CC ! KT 5-MAR-95 ! CC ! ! CC cc TP, 17-Dec-1996: geaendert fuer Alpha Prozessoren; neue SPLINE DCSPLN cc Routine aus der CERNLIB (alte nicht mehr verfuegbar) cc cc TP, 21-Nov-2001: use ranlux form Mathlib (Cernlib) as random number cc generator. cc CC C%CREATE DEDX NOKEYS C--------------------------------------------------------------- C C REAL FUNCTION DEDX(ENERG) C C THIS FUNCTION CALCULATES THE STOPPING POWER FOR C C A GIVEN ELEMENT WITH Z = NZ C C ENERG ENERGY OF THE INCIDENT PARICLE IN MEV C C NZ Z OF THE TARGET ELEMENT C C DEDX STOPPING POWER IN EV/10**15ATOMS/CM**2 C C SMASS MASS OF THE INCIDENT PARICLE IN MEV/C**2 C C SMASS, NZ ARE IN THE COMMON BLOCK DEDXPA C C COMMON /DEDXPA/ SMASS, NZ C C REFERENCE: HYDROGEN STOPPING POWERS AND RANGES IN ALL C ELEMENTS,H.H.ANDERSEN,J.F.ZIEGLER C VOLUME 3 OF THE STOPPING AND RANGES OF IONS C IN MATTER;PERGAMON PRESS 1977 C C 20-MAY-81 M.W.GLADISCH C C--------------------------------------------------------------- REAL FUNCTION DEDX(ENERG) DIMENSION A(12,92) REAL H (12),HE(12),LI(12),BE(12),B (12), 1 C (12),N (12),O (12),F (12),NE(12),NA(12),MG(12), 2 AL(12),SI(12),P (12),S (12),CL(12),AR(12),K (12), 3 CA(12),SC(12),TI(12),V (12),CR(12),MN(12),FE(12), 4 CO(12),NI(12),CU(12),ZN(12),GA(12),GE(12),AS(12), 5 SE(12),BR(12),KR(12),RB(12),SR(12),Y (12),ZR(12), 6 NB(12),MO(12),TC(12),RU(12),RH(12),PD(12),AG(12), 7 CD(12),IN(12),SN(12),SB(12),TE(12),I (12),XE(12), 8 CS(12),BA(12),LA(12),CE(12),PR(12),ND(12),PM(12), 9 SM(12),EU(12),GD(12),TB(12),DY(12),HO(12),ER(12), 1 TM(12),YB(12),LU(12),HF(12),TA(12),W (12),RE(12), 2 OS(12),IR(12),PT(12),AU(12),HG(12),TL(12),PB(12), 3 BI(12),PO(12),AT(12),RN(12),FR(12),RA(12),AC(12), 4 TH(12),PA(12),U (12) EQUIVALENCE 1 (A(1, 1),H (1)),(A(1, 2),HE(1)),(A(1, 3),LI(1)), 2 (A(1, 4),BE(1)),(A(1, 5),B (1)),(A(1, 6),C (1)), 3 (A(1, 7),N (1)), 4 (A(1, 8),O (1)),(A(1, 9),F (1)),(A(1,10),NE(1)), 5 (A(1,11),NA(1)), 6 (A(1,12),MG(1)),(A(1,13),AL(1)),(A(1,14),SI(1)), 7 (A(1,15),P (1)), 8 (A(1,16),S (1)),(A(1,17),CL(1)),(A(1,18),AR(1)), 9 (A(1,19),K (1)), 1 (A(1,20),CA(1)),(A(1,21),SC(1)),(A(1,22),TI(1)), 2 (A(1,23),V (1)), 3 (A(1,24),CR(1)),(A(1,25),MN(1)),(A(1,26),FE(1)), 4 (A(1,27),CO(1)), 5 (A(1,28),NI(1)),(A(1,29),CU(1)),(A(1,30),ZN(1)), 6 (A(1,31),GA(1)), 7 (A(1,32),GE(1)),(A(1,33),AS(1)),(A(1,34),SE(1)), 8 (A(1,35),BR(1)), 9 (A(1,36),KR(1)),(A(1,37),RB(1)),(A(1,38),SR(1)) EQUIVALENCE 1 (A(1,39),Y (1)), 2 (A(1,40),ZR(1)),(A(1,41),NB(1)),(A(1,42),MO(1)), 3 (A(1,43),TC(1)), 4 (A(1,44),RU(1)),(A(1,45),RH(1)),(A(1,46),PD(1)), 5 (A(1,47),AG(1)), 6 (A(1,48),CD(1)),(A(1,49),IN(1)),(A(1,50),SN(1)), 7 (A(1,51),SB(1)), 8 (A(1,52),TE(1)),(A(1,53),I (1)),(A(1,54),XE(1)), 9 (A(1,55),CS(1)), 1 (A(1,56),BA(1)),(A(1,57),LA(1)),(A(1,58),CE(1)), 2 (A(1,59),PR(1)), 3 (A(1,60),ND(1)),(A(1,61),PM(1)),(A(1,62),SM(1)), 4 (A(1,63),EU(1)), 5 (A(1,64),GD(1)),(A(1,65),TB(1)),(A(1,66),DY(1)), 6 (A(1,67),HO(1)), 7 (A(1,68),ER(1)),(A(1,69),TM(1)),(A(1,70),YB(1)), 8 (A(1,71),LU(1)) EQUIVALENCE 1 (A(1,72),HF(1)),(A(1,73),TA(1)),(A(1,74),W (1)), 2 (A(1,75),RE(1)), 3 (A(1,76),OS(1)),(A(1,77),IR(1)),(A(1,78),PT(1)), 4 (A(1,79),AU(1)), 5 (A(1,80),HG(1)),(A(1,81),TL(1)),(A(1,82),PB(1)), 6 (A(1,83),BI(1)), 7 (A(1,84),PO(1)),(A(1,85),AT(1)),(A(1,86),RN(1)), 8 (A(1,87),FR(1)), 9 (A(1,88),RA(1)),(A(1,89),AC(1)),(A(1,90),TH(1)), 1 (A(1,91),PA(1)), 2 (A(1,92),U (1)) COMMON / DEDXPA / SMASS, NZ, NEL, WMOL, NZL(10), ATOM(10) c COMMON / DEDXPA / SMASS, NZ DATA H /1.262,1.44,242.6,1.2E4,.1159,.0005099,5.436E4, 1 -5.052,2.049,-.3044,.01966,-.0004659/ DATA HE /1.229,1.397,484.5,5873.,.05225,.00102,2.451E4, 1 -2.158,.8278,-.1172,.007259,-.000166/ DATA LI /1.411,1.6,725.6,3013.,.04578,.00153,2.147E4, 1 -.5831,.562,-.1183,.009298,-.0002498/ DATA BE /2.248,2.59,966.,153.8,.03475,.002039,1.63E4, 1 .2779,.1745,-.05684,.005155,-.0001488/ DATA B /2.474,2.815,1206.,1060.,.02855,.002549,1.345E4, 1 -2.445,1.283,-.2205,.0156,-.000393/ DATA C /2.631,2.989,1445.,957.2,.02819,.003059,1.322E4, 1 -4.38,2.044,-.3283,.02221,-.0005417/ DATA N /2.954,3.35,1683.,1900.,.02513,.003569,1.179E4, 1 -5.054,2.325,-.3713,.02506,-.0006109/ DATA O /2.652,3.,1920.,2000.,.0223,.004079,1.046E4, 1 -6.734,3.019,-.4748,.03171,-.0007669/ DATA F /2.085,2.352,2157.,2634.,.01816,.004589,8517., 1 -5.571,2.449,-.3781,.02483,-.0005919/ DATA NE /1.951,2.199,2393.,2699.,.01568,.005099,7353., 1 -4.408,1.879,-.2814,.01796,-.0004168/ DATA NA /2.542,2.869,2628.,1854.,.01472,.005609,6905., 1 -4.959,2.073,-.3054,.01921,-.0004403/ DATA MG /3.792,4.293,2862.,1009.,.01397,.006118,6551., 1 -5.51,2.266,-.3295,.02047,-.0004637/ DATA AL /4.154,4.739,2766.,164.5,.02023,.006628,6309., 1 -6.061,2.46,-.3535,.02173,-.0004871/ DATA SI /4.15,4.7,3329.,550.,.01321,.007138,6194., 1 -6.294,2.538,-.3628,.0222,-.0004956/ DATA P /3.232,3.647,3561.,1560.,.01267,.007648,5942., 1 -6.527,2.616,-.3721,.02267,-.000504/ DATA S /3.447,3.891,3792.,1219.,.01211,.008158,5678., 1 -6.761,2.694,-.3814,.02314,-.0005125/ DATA CL /5.047,5.714,4023.,878.6,.01178,.008668,5524., 1 -6.994,2.773,-.3907,.02361,-.0005209/ DATA AR /5.731,6.5,4253.,530.,.01123,.009178,5268., 1 -7.227,2.851,-.4,.02407,-.0005294/ DATA K /5.151,5.833,4482.,545.7,.01129,.009687,5295., 1 -7.44,2.923,-.4094,.02462,-.0005411/ DATA CA /5.521,6.252,4710.,553.3,.01112,.0102,5214., 1 -7.653,2.995,-.4187,.02516,-.0005529/ DATA SC /5.201,5.884,4938.,560.9,.009995,.01071,4688., 1 -8.012,3.123,-.435,.02605,-.0005707/ DATA TI /4.862,5.496,5165.,568.5,.009474,.01122,4443., 1 -8.371,3.251,-.4513,.02694,-.0005886/ DATA V /4.48,5.055,5391.,952.3,.009117,.01173,4276., 1 -8.731,3.379,-.4676,.02783,-.0006064/ DATA CR /3.983,4.489,5616.,1336.,.008413,.01224,3946., 1 -9.09,3.507,-.4838,.02872,-.0006243/ DATA MN /3.469,3.907,5725.,1461.,.008829,.01275,3785., 1 -9.449,3.635,-.5001,.02961,-.0006421/ DATA FE /3.519,3.963,6065.,1243.,.007782,.01326,3650., 1 -9.809,3.763,-.5164,.0305,-.00066/ DATA CO /3.14,3.535,6288.,1372.,.007361,.01377,3453., 1 -10.17,3.891,-.5327,.03139,-.0006779/ DATA NI /3.553,4.004,6205.,555.1,.008763,.01428,3297., 1 -10.53,4.019,-.549,.03229,-.0006957/ DATA CU /3.696,4.175,4673.,387.8,.02188,.01479,3174., 1 -11.18,4.252,-.5791,.03399,-.0007314/ DATA ZN /4.21,4.75,6953.,295.2,.006809,.0153,3194., 1 -11.57,4.394,-.598,.03506,-.0007537/ DATA GA /5.041,5.697,7173.,202.6,.006725,.01581,3154., 1 -11.95,4.537,-.6169,.03613,-.0007759/ DATA GE /5.554,6.3,6496.,110.,.009689,.01632,3097., 1 -12.34,4.68,-.6358,.03721,-.0007981/ DATA AS /5.323,6.012,7611.,292.5,.006447,.01683,3024., 1 -12.72,4.823,-.6547,.03828,-.0008203/ DATA SE /5.874,6.656,7395.,117.5,.007684,.01734,3006., 1 -13.11,4.965,-.6735,.03935,-.0008425/ DATA BR /5.611,6.335,8046.,365.2,.006244,.01785,2928., 1 -13.4,5.083,-.6906,.04042,-.0008675/ DATA KR /6.411,7.25,8262.,220.,.006087,.01836,2855., 1 -13.69,5.2,-.7076,.0415,-.0008925/ DATA RB /5.694,6.429,8478.,292.9,.006087,.01886,2855., 1 -13.92,5.266,-.714,.04173,-.0008943/ DATA SR /6.339,7.159,8693.,330.3,.006003,.01937,2815., 1 -14.14,5.331,-.7205,.04196,-.0008962/ DATA Y /6.407,7.234,8907.,367.8,.005889,.01988,2762., 1 -14.36,5.397,-.7269,.04219,-.000898/ DATA ZR /6.734,7.603,9120.,405.2,.005765,.02039,2704., 1 -14.59,5.463,-.7333,.04242,-.0008998/ DATA NB /6.902,7.791,9333.,442.7,.005587,.0209,2621., 1 -16.22,6.094,-.8225,.04791,-.001024/ DATA MO /6.425,7.248,9545.,480.2,.005367,.02141,2517., 1 -17.85,6.725,-.9116,.05339,-.001148/ DATA TC /6.799,7.671,9756.,517.6,.005315,.02192,2493., 1 -17.96,6.752,-.9135,.05341,-.001147/ DATA RU /6.108,6.887,9966.,555.1,.005151,.02243,2416., 1 -18.07,6.779,-.9154,.05342,-.001145/ DATA RH /5.924,6.677,1.018E4,592.5,.004919,.02294,2307., 1 -18.18,6.806,-.9173,.05343,-.001143/ DATA PD /5.238,5.9,1.038E4,630.,.004758,.02345,2231., 1 -18.28,6.833,-.9192,.05345,-.001142/ DATA AG /5.623,6.354,7160.,337.6,.01394,.02396,2193., 1 -18.39,6.86,-.9211,.05346,-.00114/ DATA CD /5.814,6.554,1.08E4,355.5,.004626,.02447,2170., 1 -18.62,6.915,-.9243,.0534,-.001134/ DATA IN /6.23,7.024,1.101E4,370.9,.00454,.02498,2129., 1 -18.85,6.969,-.9275,.05335,-.001127/ DATA SN /6.41,7.227,1.121E4,386.4,.004474,.02549,2099., 1 -19.07,7.024,-.9308,.05329,-.001121/ DATA SB /7.5,8.48,8608.,348.,.009074,.026,2069., 1 -19.57,7.225,-.9603,.05518,-.001165/ DATA TE /6.979,7.871,1.162E4,392.4,.004402,.02651,2065., 1 -20.07,7.426,-.9899,.05707,-.001209/ DATA I /7.725,8.716,1.183E4,394.8,.004376,.02702,2052., 1 -20.56,7.627,-1.019,.05896,-.001254/ DATA XE /8.231,9.289,1.203E4,397.3,.004384,.02753,2056., 1 -21.06,7.828,-1.049,.06085,-.001298/ DATA CS /7.287,8.218,1.223E4,399.7,.004447,.02804,2086., 1 -20.4,7.54,-1.004,.05782,-.001224/ DATA BA /7.899,8.911,1.243E4,402.1,.004511,.02855,2116., 1 -19.74,7.252,-.9588,.05479,-.001151/ DATA LA /8.041,9.071,1.263E4,404.5,.00454,.02906,2129., 1 -19.08,6.964,-.9136,.05176,-.001077/ DATA CE /7.489,8.444,1.283E4,406.9,.00442,.02957,2073., 1 -18.43,6.677,-.8684,.04872,-.001003/ DATA PR /7.291,8.219,1.303E4,409.3,.004298,.03008,2016., 1 -17.77,6.389,-.8233,.04569,-.0009292/ DATA ND /7.098,8.,1.323E4,411.8,.004182,.03059,1962., 1 -17.11,6.101,-.7781,.04266,-.0008553/ DATA PM /6.91,7.786,1.343E4,414.2,.004058,.0311,1903., 1 -16.45,5.813,-.733,.03963,-.0007815/ DATA SM /6.728,7.58,1.362E4,416.6,.003976,.03161,1865., 1 -15.79,5.526,-.6878,.0366,-.0007077/ DATA EU /6.551,7.38,1.382E4,419.,.003877,.03212,1819., 1 -15.13,5.238,-.6426,.03357,-.0006339/ DATA GD /6.739,7.592,1.402E4,421.4,.003863,.03263,1812., 1 -14.47,4.95,-.5975,.03053,-.0005601/ DATA TB /6.212,6.996,1.421E4,423.9,.003725,.03314,1747., 1 -14.56,4.984,-.6022,.03082,-.0005668/ DATA DY /5.517,6.21,1.44E4,426.3,.003632,.03365,1703., 1 -14.65,5.018,-.6069,.03111,-.0005734/ DATA HO /5.219,5.874,1.46E4,428.7,.003498,.03416,1640., 1 -14.74,5.051,-.6117,.03141,-.0005801/ DATA ER /5.071,5.706,1.479E4,433.,.003405,.03467,1597., 1 -14.83,5.085,-.6164,.0317,-.0005867/ DATA TM /4.926,5.542,1.498E4,433.5,.003342,.03518,1567., 1 -14.91,5.119,-.6211,.03199,-.0005933/ DATA YB /4.787,5.386,1.517E4,435.9,.003292,.03569,1544., 1 -15.,5.153,-.6258,.03228,-.0006/ DATA LU /4.893,5.505,1.536E4,438.4,.003243,.0362,1521., 1 -15.09,5.186,-.6305,.03257,-.0006066/ DATA HF /5.028,5.657,1.555E4,440.8,.003195,.03671,1499., 1 -15.18,5.22,-.6353,.03286,-.0006133/ DATA TA /4.738,5.329,1.574E4,443.2,.003186,.03722,1494., 1 -15.27,5.254,-.64,.03315,-.0006199/ DATA W /4.574,5.144,1.593E4,442.4,.003144,.03773,1475., 1 -15.67,5.392,-.6577,.03418,-.0006426/ DATA RE /5.2,5.851,1.612E4,441.6,.003122,.03824,1464., 1 -16.07,5.529,-.6755,.03521,-.0006654/ DATA OS /5.07,5.704,1.63E4,440.9,.003082,.03875,1446., 1 -16.47,5.667,-.6932,.03624,-.0006881/ DATA IR /4.945,5.563,1.649E4,440.1,.002965,.03926,1390., 1 -16.88,5.804,-.711,.03727,-.0007109/ DATA PT /4.476,5.034,1.667E4,439.3,.002871,.03977,1347., 1 -17.28,5.942,-.7287,.0383,-.0007336/ DATA AU /4.856,5.46,1.832E4,438.5,.002542,.04028,1354., 1 -17.02,5.846,-.7149,.0374,-.0007114/ DATA HG /4.308,4.843,1.704E4,487.8,.002882,.04079,1352., 1 -17.84,6.183,-.7659,.04076,-.0007925/ DATA TL /4.723,5.311,1.722E4,537.,.002913,.0413,1366., 1 -18.66,6.52,-.8169,.04411,-.0008737/ DATA PB /5.319,5.982,1.74E4,586.3,.002871,.04181,1347., 1 -19.48,6.857,-.8678,.04747,-.0009548/ DATA BI /5.956,6.7,1.78E4,677.,.00266,.04232,1336., 1 -19.55,6.871,-.8686,.04748,-.0009544/ DATA PO /6.158,6.928,1.777E4,586.3,.002812,.04283,1319., 1 -19.62,6.884,-.8694,.04748,-.000954/ DATA AT /6.204,6.979,1.795E4,586.3,.002776,.04334,1302., 1 -19.69,6.898,-.8702,.04749,-.0009536/ DATA RN /6.181,6.954,1.812E4,586.3,.002748,.04385,1289., 1 -19.76,6.912,-.871,.04749,-.0009532/ DATA FR /6.949,7.82,1.83E4,586.3,.002737,.04436,1284., 1 -19.83,6.926,-.8718,.0475,-.0009528/ DATA RA /7.506,8.448,1.848E4,586.3,.002727,.04487,1279., 1 -19.9,6.94,-.8726,.04751,-.0009524/ DATA AC /7.649,8.609,1.866E4,586.3,.002697,.04538,1265., 1 -19.97,6.953,-.8733,.04751,-.000952/ DATA TH /7.71,8.679,1.883E4,586.3,.002641,.04589,1239., 1 -20.04,6.967,-.8741,.04752,-.0009516/ DATA PA /7.407,8.336,1.901E4,586.3,.002603,.0464,1221., 1 -20.11,6.981,-.8749,.04752,-.0009512/ DATA U /7.29,8.204,1.918E4,586.3,.002573,.04691,1207., 1 -20.18,6.995,-.8757,.04753,-.0009508/ C ENER=ENERG*931501.6/SMASS IF ( ENER . LT . 10. ) GO TO 100 IF ( ENER . LT . 1000.) GO TO 200 C C REDUCED ENERGY LARGER THAN 1 MEV C R = 0. DO 10 III=1,5 R = R + A(III+7,NZ)*ALOG(ENER)**FLOAT(III-1) 10 CONTINUE VEL2=(ENERG**2+2.*ENERG*SMASS)/(ENERG+SMASS)**2 c-lz DEDX=A(6,NZ)/VEL2*ALOG(A(7,NZ)*VEL2/(1.-VEL2)-VEL2-R)!Original falsch DEDX=A(6,NZ)/VEL2*(ALOG(A(7,NZ)*VEL2/(1.-VEL2))-VEL2-R) GO TO 500 C C REDUCED ENERGY LESS THAN 10. KEV C 100 DEDX=A(1,NZ)*SQRT(ENER) GO TO 500 C C REDUCED ENERGY BETWEEN 10 KEV AND 1 MEV C 200 DEDX=1./(1./(A(2,NZ)*ENER**.45)+1./(A(3,NZ)/ENER* 1 ALOG(1.+A(4,NZ)/ENER+A(5,NZ)*ENER))) 500 CONTINUE RETURN END C%CREATE DEDXM NOKEYS C----------------------------------------------------------------------- C C NAME: DEDXM C C FUNKTION: BERECHNET DIE DEDX-FUNKTION FUER BELIEBIGE GEMISCHE C C AUFRUF: DEDXV = DEDXM ( EKIN ) C C PARAMETER: C C EKIN R*4 KIN. ENERGIE DES TEILCHENS IN MEV C C FUNKTIONSWERT: C C DEDXV R*4 ENERGIEVERLUST IN MEV/(G/CM**2)) C C BEMERKUNG: DIE DATEN MUESSEN AUF DEM COMMON-BLOCK /DEDXPA/ C WIE FOLGT ABGELEGT WORDEN SEIN: C C SMASS R*4 MASSE DER TEILCHENS IN MEV/C**2 C NZ I*4 KERNLADUNGSZAHL FUER "DEDX"-FUNKTION C NEL I*4 ZAHL DER ELEMENTE C WMOL R*4 "MOLEKULARGEWICHT" = SUMME ( ATOM(I)*AWT(I) ) C NZL R*4 KERNLADUNGSZAHLEN DER ELEMENTE C ATOM R*4 ZAHL DER ATOME DER ELEMENTE C C----------------------------------------------------------------------- REAL FUNCTION DEDXM (EKIN) EXTERNAL DEDX REAL SDEDX, EKIN COMMON / DEDXPA / SMASS, NZ, NEL, WMOL, NZL(10), ATOM(10) C---- SDEDX = 0.0 DO 10 I=1,NEL NZ = NZL(I) SDEDX = SDEDX + ( DEDX(EKIN) * ATOM(I) ) 10 CONTINUE DEDXM = SDEDX * 602.2045 / WMOL RETURN END C%CREATE DGESMP NOKEYS C----------------------------------------------------------------- C C NAME: DGESMP C C FUNKTION: ERSTELLT DIE REICHWEITE-ENERGIE TABELLE C DURCH INTEGRATION DER DX/DE-FUNKTION. C C AUFRUF: CALL DGESMP(F,A,B,ACC,ANS,ERROR,IFLAG) C C PARAMETER: C C F R*4 FUNCTION MIT EINEM ARGUMENT c c F R*8 besser wegen Alpha Prozessor; auch Argument muss c als REAL*8 uebergeben werden, sonst gibt's c Probleme bei der Wertzuweisung; TP, 17-Dec-1996 c C A R*4 START-WERT C B R*4 END-WERT C ACC R*4 RELATIVE GENAUIGKEIT C ANS R*4 WERT DES INTEGRALS C ERROR R*4 GESCHAETZTER FEHLER C IFLAG I*4 FEHLER-PARAMETER C =1 O.K. C =2 ZU VIELE UNTERINTERVALLE C =3 ZU VIELE FUNKTIONSAUFRUFE C =4 MEHR ALS 200 ZWISCHENPUNKTE C C BEMERKUNG: DIESE PROGRAMM UNTERSCHEIDET SICH IN FOLGENDEN C PUNKTEN VOM BISHERIGEN "SIMP"-PROGRAMM: C C (1) DER RELATIVE FEHLER WIRD BEI JEDER HALBIERUNG C DES INTERVALLS EBENFALLS HALBIERT UND NICHT C WIE ZUVOR DURCH 1.4 DIVIDIERT ! C C (2) DIE INTEGRATION SELBST WURDE ZUR SICHERHEIT C DOPPELTGENAU VORGENOMMEN. DIE FUNKTION "F" MUSS C DOPPELTGENAU BERECHNET WERDEN. C C (3) DIE TABELLE REICHWEITE ALS FUNKTION DER ENERGIE C WIRD ERSTELLT ( AUF /WORK/ ). C ======>> ALS REAL*8 WERTE C C-------------------------------------------------------------------- SUBROUTINE DGESMP(F,A,B,ACC,ANS,ERROR,IFLAG) IMPLICIT REAL*8 (A-H,O-Z) real*8 f EXTERNAL F REAL*4 A,B,ACC,ANS,ERROR,AREA REAL*8 XVAL,YVAL,COEFF c-kt Deklaration neuer Bezeichner für Ein- und Ausgabeunit: SYSOUT, SYSIN INTEGER*2 SYSOUT, SYSIN c-kt Deklaration der logischen Variable für Tests LOGICAL TEST DIMENSION FV(5),LORR(50),F1T(50),F2T(50),F3T(50),DAT(50), 1 ARESTT(50),ESTT(50),EPST(50),PSUM(50), 2 XL(50),XR(50) c-kt Deklaration neuer Bezeichner für Ein- und Ausgabeunit: SYSOUT, SYSIN COMMON 1 / STATUS / SYSOUT, SYSIN, TEST COMMON / WORK / NVMAX,NVAL,XVAL(200),YVAL(200),COEFF(200,3) c COMMON / WORK / NVMAX,NVAL,XVAL(200),YVAL(200) DATA NTMAX/50/ C---- NVAL = 1 XVAL(NVAL) = A YVAL(NVAL) = 0. C---- U = 9.0D-10 FOURU = 4.0D00 * U IFLAG = 1 c-kt Hier wird die relative Genauigkeit der Interpolation zugewiesen: EPS = ACC ERROR = 0.0 LVL = 1 LORR(LVL) = 1 PSUM(LVL) = 0.0D00 sum = 0.0d00 C---- IF( TEST ) THEN WRITE (SYSOUT,1802) LVL,LORR(LVL),PSUM(LVL),SUM,XL(LVL), 1 XR(LVL),ARESTT(LVL),ESTT(LVL),EPST(LVL) 1802 FORMAT(' ',' # T',2X,' SUMME ',2X,' TOTAL ', 1 2X,' LINKS ', 2 2X,' RECHTS ',2X,'|F|-INTEGRAL', 3 2X,' F-INTEGRAL ',2X,' EPSILON '// 4 ' ',I2,1X,I1,2X,7(E12.5,2X)) ENDIF C---- ALPHA = A DA = B - A AREA = 0.0 AREST = 0.0 FV(1) = F(ALPHA) FV(3) = F(ALPHA+0.5D00*DA) FV(5) = F(ALPHA+DA) KOUNT = 3 WT = DA / 6.0D00 EST = WT * ( FV(1) + 4.0D00 * FV(3) + FV(5)) C C=================== STARTE ALGORITHMUS C 1 DX = 0.5D00*DA FV(2) = F(ALPHA+0.5D00*DX) FV(4) = F(ALPHA+1.5D00*DX) KOUNT = KOUNT + 2 WT = DX / 6.0D00 ESTL = WT * ( FV(1) + 4.0D00 * FV(2) + FV(3) ) ESTR = WT * ( FV(3) + 4.0D00 * FV(4) + FV(5) ) SUM = ESTL + ESTR ARESTL = WT * ( ABS(FV(1)) + ABS(FV(2)*4.0D00) + ABS(FV(3))) ARESTR = WT * ( ABS(FV(3)) + ABS(FV(4)*4.0D00) + ABS(FV(5))) AREA = AREA + ((ARESTL+ARESTR)-AREST) DIFF = EST - SUM C C=================== GENAUIGKEIT ERREICHT ?? C c-kt Diese Bedingung ist für z.B. Aluminium viel zu früh erfüllt: IF ( ABS(DIFF) .LE. EPS*ABS(AREA) ) GOTO 2 C C=================== RECHENGENAUIGKEIT ERREICHT ?? C IF ( ABS(DX) .LE. FOURU*ABS(ALPHA) ) THEN IFLAG = 2 ELSE C C=================== ZU VIELE UNTERINTERVALLE C IF ( LVL .GE. NTMAX ) THEN IFLAG = 2 ELSE C C=================== ZU VIELE FUNKTIONS-AUFRUFE ?? C IF ( KOUNT .GE. 2000 ) THEN IFLAG = 3 ELSE GOTO 11 ENDIF ENDIF ENDIF C=================== BRECHE DIE INTEGRATION DES ABSCHNITTES AB ! GOTO 2 C C=================== INTEGRAL NOCH NICHT GENAU GENUG C RETTE DIE RECHTE HAELFTE AUF (LVL) C 11 LVL = LVL + 1 LORR(LVL) = 0 PSUM(LVL) = 0.D00 F1T(LVL) = FV(3) F2T(LVL) = FV(4) F3T(LVL) = FV(5) XL (LVL) = ALPHA + DX XR (LVL) = XL (LVL) + DX DA = DX DAT(LVL) = DX AREST = ARESTL ARESTT(LVL) = ARESTR EST = ESTL ESTT(LVL) = ESTR EPS = EPS / 2.D00 EPST(LVL) = EPS FV(5) = FV(3) FV(3) = FV(2) C---- IF ( TEST ) THEN WRITE (SYSOUT,1801) LVL,LORR(LVL),PSUM(LVL),SUM,XL(LVL), 1 XR(LVL),ARESTT(LVL),ESTT(LVL),EPST(LVL) 1801 FORMAT(' ',I2,1X,I1,2X,7(E12.5,2X)) ENDIF C---- GOTO 1 C C=================== GENAUIGKEIT ERREICHT C 2 ERROR = ERROR + DIFF/15.0D00 3 IF ( LORR(LVL) .EQ. 0 ) GOTO 4 C C=================== VORHERIGES INTERVALL SCHON BERECHNET C SUM = PSUM(LVL) + SUM IF ( TEST ) WRITE (SYSOUT,1801) LVL,LORR(LVL),PSUM(LVL),SUM, 1 XL(LVL),XR(LVL),ARESTT(LVL), 2 ESTT(LVL),EPST(LVL) LVL = LVL - 1 IF ( LVL .GT. 1 ) GOTO 3 C C=================== FERTIG C ANS = SUM NVAL = NVAL + 1 XVAL(NVAL) = B YVAL(NVAL) = ANS IF ( TEST ) THEN WRITE (SYSOUT,1701) (KK,XVAL(KK),YVAL(KK),KK=1,NVAL) 1701 FORMAT(' ',' N = ',I3,' X = ',E12.5,' Y = ',E12.5) WRITE (SYSOUT,1900) SUM,IFLAG,ERROR 1900 FORMAT(/' INTEGRAL=',E12.5,2X,'IFLAG=',I1,2X,'ERROR=',E12.5/) ENDIF RETURN C C=================== NAECHSTES INTERVALL ZU BERECHNEN C 4 PSUM(LVL) = SUM LORR(LVL) = 1 IF ( NVAL .GE. NVMAX ) GOTO 999 C---- ****** EXKURS: ERSTELLE TEILINTEGRALE NVAL = NVAL + 1 XVAL(NVAL) = ALPHA + DA DUMMY = 0.0D00 DO 22 KK=1,LVL IF ( LORR(KK) .EQ. 1 ) DUMMY = DUMMY + PSUM(KK) 22 CONTINUE YVAL(NVAL) = DUMMY C---- ALPHA = ALPHA + DA DA = DAT(LVL) FV(1) = F1T(LVL) FV(3) = F2T(LVL) FV(5) = F3T(LVL) AREST = ARESTT(LVL) EST = ESTT(LVL) EPS = EPST(LVL) C---- IF ( TEST ) WRITE (SYSOUT,1801) LVL,LORR(LVL),PSUM(LVL),SUM, 1 XL(LVL),XR(LVL),ARESTT(LVL), 2 ESTT(LVL),EPST(LVL) C---- GOTO 1 C MEHR ALS 200 PUNKTE 999 IFLAG = 4 RETURN END C%CREATE DGETAB NOKEYS C----------------------------------------------------------------------- C C NAME: DGETAB C C FUNKTION: BERECHNET EINE REICHWEITE-ENERGIE TABELLE UND C INTERPOLIERT DIESE MIT EINER KUBISCHEN SPLINE- C FUNKTION. DIE KOEFFIZIENTEN WERDEN AUF DEM FELD C RGDATA ABGELEGT. C MIT "RGEVAL" KANN ANSCHLIESSEND BEI GEGEBENER C SCHICHTDICKE UND ANFANGSENERGIE DIE ENDENERGIE C DES TEILCHENS BERECHNET WERDEN. C C AUFRUF: CALL DGETAB ( PMASS,E1,E2,PREC, C NZ,AWT,NATOM,NE, C RGDATA,NRG, C IFLAG ) C C PARAMETER: C C PMASS R*4 MASSE DES TEILCHENS IN MEV/C2 C E1 R*4 MINDESTENERGIE IN MEV C E2 R*4 MAXIMALENERGIE IN MEV C PREC R*4 RELATIVE GENAUIGKEIT DER INTERPOLATION C C NZ(NE) I*4 KERNLADUNGSZAHL C AWT(NE) R*4 ATOMGEWICHT IN GRAMM C NATOM(NE) I*4 ZAHL DER ATOME IM MOLEKUEL C NE I*4 ZAHL DER ELEMENTE C C RGDATA(5,NRG) R*8 SPLINE COEFFIZIENTEN C NRG I*4 ANZAHL DER PUNKTE C C IFLAG I*4 FEHLERNUMMER; =0 KEIN FEHLER C C BEMERKUNG: DAS PROGRAMM BENUTZT EINEN COMMON-BLOCK /WORK/ C ALS ARBEITSSPEICHER. C C----------------------------------------------------------------------- SUBROUTINE DGETAB ( PMASS, E1, E2, PREC, 1 NZ, AWT, NATOM, NE, 2 RGDATA, NRG, 3 IFLAG ) C---- c real*8 indedx EXTERNAL INDEDX c INTEGER NZ(NE), NE, NRG, IFLAG, NATOM(NE) c c-kt Deklaration neuer Bezeichner für Ein- und Ausgabeunit: SYSOUT, SYSIN c INTEGER*4 TEST INTEGER*2 SYSOUT, SYSIN REAL PMASS, PREC, AWT(NE) c REAL*4 X(200),Y(200),DER(200,2),Z(200),FVAL(200),FDER(200,2) c INTEGER IOP/-1/,MSP/0/ REAL*8 RGDATA(5,NRG), XVAL,YVAL,COEFF c c-kt Deklaration eines neuen Bezeichners für die Ausgabeunit: SYSOUT c COMMON 1 / STATUS / SYSOUT, SYSIN, TEST COMMON / WORK / NVMAX,NVAL,XVAL(200),YVAL(200),COEFF(200,3) COMMON / DEDXPA / SMASS, NZLL, NEL, WMOL, NZL(10), ATOM(10) COMMON / SPAPPR / SECD1,SECDN,VOFINT,IERR,NXY C---- c c for new spline routine DCSPLN c parameter (M = 1) integer*4 nsub c real*8 x(200), y(200,m), a(200,m), b(200,m), c(200,m), d(200,m) c c--------------------------------------------------------------------------- c NVMAX = 200 C---- C UEBERGEBE TEILCHEN- UND MATERIAL-PARAMETER AN INDEDX C---- SMASS = PMASS NEL = NE WMOL = 0.0 DO 10 I=1,NE NZL(I) = NZ(I) WMOL = WMOL + ( FLOAT(NATOM(I)) * AWT(I) ) ATOM(I) = FLOAT(NATOM(I)) 10 CONTINUE C---- C ERMITTLE DIE STUETZPUNKTE DES SPLINE-POLYNOMS C---- EKMIN = E1 EKMAX = E2 PREC0 = PREC DO 15 NCOUNT=1,20 CALL DGESMP(INDEDX,EKMIN,EKMAX,PREC0,ANS,ERROR,IFLAG) IF ( IFLAG .NE. 1 .OR. NVAL .GT. NRG ) THEN c-lz it seems that it lacks something hier,but what? 20/12/91 L.Zhang c-kt Ich gehe nicht mehr davon aus, daß hier etwas fehlt. KT 10-JAN-95 ELSE IF ( NVAL .LE. 40 ) GOTO 18 ENDIF c-kt Wenn 20*Verringern der rel. Genauigkeit nichts bringt, ... PREC0 = PREC0 * 1.414 15 CONTINUE c-kt ... wird IFLAG=5 gesetzt, ... IFLAG = 5 c-kt ... was wir hier dem Benutzer mitteilen wollen (Fortsetzung in MCV3K): WRITE (SYSOUT,5555) PREC0 5555 FORMAT(' Der Parameter von ''BUILD_TAB'' sollte auf ',E7.1, 1 ' gesetzt werden für die Schicht') GOTO 100 C---- C DRUCKE DIE TABELLE AUS C---- 18 WRITE (SYSOUT,1702) 1702 FORMAT(' ENERGIE-REICHWEITE TABELLE'/) WRITE (SYSOUT,1701) (KK,XVAL(KK),YVAL(KK),KK=1,NVAL) 1701 FORMAT(' ',' N = ',I3,' X = ',1PE12.5,' Y = ',1PE12.5) C WRITE (SYSOUT,1900) ANS,IFLAG,PREC0 1900 FORMAT(/' INTEGRAL=',1PE12.5,2X,'IFLAG=',I1,2X,'ERROR=',1PE12.5/) C---- C ERSTELLE DIE SPLINE-KOEFFIZIENTEN C---- DO ISP = 1,NVAL X(ISP) = XVAL(ISP) Y(ISP,m) = YVAL(ISP) END DO c c new routine double precision DCSPLN; the old one SPLIN3 is obsolete; c unfortunately other parameters are needed; c c number of subintervals must be at least 2 c c the spline is a polynomial a + b*x + c*x**2 + d*x**3 c c c 1. parameter: nsub = nval-1 : number of subintervals [x(i-1),x(i)] c 2. parameter: x : one dimensional array of at least nsub+1 = nval c elements; contains abscissa values. Unchanged on c return c 3. parameter: m : number of components of the vector valued spline c function; here m = 1 c 4. parameter: y : two dimensional array of dimension (0:NDIM,d) c with d>=m, here d=1; c c 5. parameter: ndim : at least nsub c c 6. parameter: mode : = 0 if natural spline f''(x0) = f''(xn) c = 1 if f''(x0)=f''(x1) and f''(xn-1) = f''(xn) c c 7. parameter : a(NDIM,D) d>=m spline coefficient zero order c 8. parameter : b(NDIM,D) d>=m spline coefficient first order c 9. parameter : c(NDIM,D) d>=m spline coefficient second order c 10. parameter: d(NDIM,D) d>=m spline coefficient third order c nsub = nval - 1 c call dcspln (nsub, x, m, y, nsub, 0, a, b, c, d) c IFLAG = IERR IF ( IFLAG .NE. 0 ) RETURN C---- C KOPIERE DIE KOEFFIZIENTEN AUF "RGDATA" C---- NRG = NVAL DO 50 L=1,NVAL RGDATA(1,L) = XVAL(L) RGDATA(2,L) = YVAL(L) IF( L .EQ. NVAL) THEN RGDATA(5,L) = RGDATA(5,L-1) ELSE RGDATA(3,L) = b(l,1) RGDATA(4,L) = c(l,1) RGDATA(5,L) = d(l,1) END IF 50 continue c do l = 1, nval c write(sysout,'(i4,1x,d10.4,1x,d10.4,1x,d10.4,1x,d10.4,1x,d10.4,1x,d10.4)') c 1 l, rgdata(1,l), rgdata(2,l), a(l,1), c 2 rgdata(3,l), rgdata(4,l), rgdata(5,l) c enddo c c c DO ISP = 1,NVAL c X(ISP) = XVAL(ISP) c Y(ISP) = YVAL(ISP) c END DO c-kt Weniger als 4 Werte in der Reichweite-Energie-Tab. verträgt SPLIN3 nicht c CALL SPLIN3(X,Y,DER,NVAL,NVMAX,Z,FVAL,FDER,MSP,NVMAX,IOP) c IFLAG = IERR c IF ( IFLAG .NE. 0 ) RETURN C---- C KOPIERE DIE KOEFFIZIENTEN AUF "RGDATA" C---- c NRG = NVAL c DO 50 L=1,NVAL c RGDATA(1,L) = XVAL(L) c RGDATA(2,L) = YVAL(L) c RGDATA(3,L) = DER(L,1) c RGDATA(4,L) = DER(L,2)/2 c IF( L .EQ. NVAL) THEN c RGDATA(5,L) = RGDATA(5,L-1) c ELSE c RGDATA(5,L) = (DER(L,2) - DER(L+1,2))/(X(L) - X(L+1))/6 c END IF c c c 50 CONTINUE C IFLAG = 0 C 100 RETURN END C%CREATE DRANGE NOKEYS C----------------------------------------------------------------------- C C NAME: DRANGE C C FUNKTION: BERECHNET DIE REICHWEITE EINES TEILCHENS MIT EINER C KINETISCHEN ENERGIE "EKIN" IN EINEM GEGEBENEN C MATERIAL. C C AUFRUF: R = DRANGE ( EKIN,RGDATA,NRG ) C C PARAMETER: C C EKIN R*4 ANFANGSENERGIE C RGDATA(5,NRG) R*8 SPLINE COEFFIZIENTEN C NRG I*4 ANZAHL DER INTERPOLATIONS-PUNKTE C C FUNKTIONSWERT: C C R R*4 REICHWEITE IN G/CM**2 C C BEM.: FALLS DIE ENERGIE GROESSER ALS DER HOECHSTE TABELLIERTE C ENERGIEWERT IST, WIRD DAS PROGRAMM MIT "STOP 12" ABGEBROCHEN. c-kt Geändert! Der Funktion wird jetzt in diesem Fall (korrekt: größer gleich) c-kt der Wert "-1." zugewiesen; die Fehlerbehandlung wird der auf- c-kt rufenden Routine überlassen. C C----------------------------------------------------------------------- REAL FUNCTION DRANGE ( EKIN, RGDATA, NRG ) REAL EKIN REAL*8 EK, DX1, RGDATA(5,NRG) INTEGER NRG, TEST c-kt Deklaration eines neuen Bezeichners für die Ausgabeunit: SYSOUT INTEGER*2 SYSOUT, SYSIN c-kt Deklaration eines neuen Bezeichners für die Ausgabeunit: SYSOUT COMMON 1 / STATUS / SYSOUT, SYSIN, TEST C---- NVAL = NRG EK = EKIN C IF ( EK .GE. RGDATA(1,NVAL) ) THEN c-kt die Fehlerbehandlung wird der aufrufenden Routine überlassen: WRITE(SYSOUT,3001) 3001 FORMAT(//' Schwerer Fehler in der Funktion ''DRANGE'':') WRITE(SYSOUT,3002) 3002 FORMAT(' ',T5,'Ein Teilchen hatte einen höheren Wert für die'// 1'kinetische Energie'/T5,'als der größte tabellierte Energiewert.') WRITE(SYSOUT,3003) EK, RGDATA(1,NVAL) 3003 FORMAT(' ',T5,'kinetische Energie des Teilchens',T50,D22.16,/T5, 1 'größter tabellierter Energiewert',T50,D22.16) DRANGE = -1. RETURN c-kt Ende c-kt WRITE(SYSOUT,1750) EK,RGDATA(1,NVAL) c-kt 1750 FORMAT(' EKIN (',E12.5,') > EMAX (',E12.5,')') c-kt STOP 12 ENDIF C---- C SUCHE START-INTERVALL C---- DO 10 K=NVAL,1,-1 IF ( EK .GT. RGDATA(1,K) ) GOTO 20 10 CONTINUE 15 DRANGE = 0.0 RETURN C---- C START-INTERVALL GEFUNDEN, BERECHNE REICHWEITE C---- 20 N1 = K DX1 = EKIN - RGDATA(1,N1) DRANGE = DX1**3 * RGDATA(5,N1) + DX1**2 * RGDATA(4,N1) + 2 DX1 * RGDATA(3,N1) + RGDATA(2,N1) RETURN END C%CREATE ELOSS2 NOKEYS C----------------------------------------------------------------------- C C NAME: ELOSS2 C C FUNKTION: BERECHNET DEN ENERGIEVERLUST EINES TEILCHENS NACH C DURCHLAUFEN EINER SCHICHT MIT DER SCHICHTDICKE C "RTHICK" IN G/CM**2 BEI EINER ANFANGSENERGIE VON C "EKIN" MEV. C C AUFRUF: DE = ELOSS2( EKIN,RTHICK,RGDATA,NRG ) C C PARAMETER: C C EKIN R*4 ANFANGSENERGIE C RTHICK R*4 SCHICHTDICKE IN G/CM**2 C RGDATA(5,NRG) R*4 SPLINE COEFFIZIENTEN C NRG I*4 ANZAHL DER INTERPOLATIONS-PUNKTE C C FUNKTIONSWERT: C C ELOSS2 R*4 ENERGIEVERLUST IN MEV C C BEM.: 1. FALLS DIE ENERGIE GROESSER ALS DER HOECHSTE TABELLIERTE C ENERGIEWERT IST, WIRD DAS PROGRAMM MIT "STOP 12" ABGEBROCHEN. c-kt Geändert! Der Funktion wird jetzt in diesem Fall (korrekt: größer gleich) c-kt der Wert "-1." zugewiesen; die Fehlerbehandlung wird der c-kt aufrufenden Routine überlassen. C 2. DAS PROGRAMM RECHNET JETZT INTERN MIT REAL*8. c-kt 3. Bislang wurde das Programm ebenfalls kommentarlos abgebrochen, c-kt wenn kein Reichweiten-Wert in der Energie-Reichweite-Tabelle c-kt kleiner ist als die nach Abzug der Schichtdicke [g/cm**2] von c-kt der Eingangs-Reichweite verbleibende Dicke, also kein Endinter- c-kt vall gefunden wurde. In diesem Fall wird jetzt der Funktion der c-kt Wert "-2." zugewiesen. c-kt 4. Auch bei einer fehlerhaften Lösung des Polynoms im Falle dreier c-kt reeller Nullstellen wurde das Programm bislang kommentarlos ab- c-kt gebrochen; jetzt wird der Funktion der Wert "-3." zugewiesen. C C----------------------------------------------------------------------- REAL FUNCTION ELOSS2 (EKIN,RTHICK,RGDATA,NRG) REAL EKIN, RTHICK REAL*8 RGDATA(5,NRG) INTEGER NRG, TEST c-kt Deklaration eines neuen Bezeichners für die Ausgabeunit: SYSOUT INTEGER*2 SYSOUT, SYSIN c-kt Deklaration eines neuen Bezeichners für die Ausgabeunit: SYSOUT COMMON 1 / STATUS / SYSOUT, SYSIN, TEST REAL*8 EK, RT, DX1, RANGE, RGRENZ, 1 X,XDUM,RHO,PHI,EKOUT,DRITTL,RDRITT,PDRITT,R,S,T,P,Q,D,U,V c-kt Die Genauigkeit der Konstante ist von 14 auf 16 Stellen erweitert: DATA DRITTL/0.3333333333333333D00/ C---- NVAL = NRG EK = EKIN RT = RTHICK C IF ( EK .GE. RGDATA(1,NVAL) ) THEN c-kt die Fehlerbehandlung wird der aufrufenden Routine überlassen: WRITE(SYSOUT,3004) 3004 FORMAT(//' Schwerer Fehler in der Funktion ''ELOSS2'':') WRITE(SYSOUT,3005) 3005 FORMAT(' ',T5,'Ein Teilchen hatte einen höheren Wert für die'// 1'kinetische Energie'/T5,'als der größte tabellierte Energiewert.') WRITE(SYSOUT,3006) EK, RGDATA(1,NVAL) 3006 FORMAT(' ',T5,'kinetische Energie des Teilchens',T50,D22.16,/T5, 1 'größter tabellierter Energiewert',T50,D22.16) ELOSS2 = -1. RETURN c-kt Ende c-kt WRITE(SYSOUT,1750) EK,RGDATA(1,NVAL) c-kt 1750 FORMAT(' EKIN (',E12.5,') > EMAX (',E12.5,')') c-kt STOP 12 ENDIF C---- C SUCHE START-INTERVALL C---- DO 10 K=NVAL,1,-1 IF ( EK .GT. RGDATA(1,K) ) GOTO 20 10 CONTINUE 15 ELOSS2 = EK RETURN C---- C START-INTERVALL GEFUNDEN, BERECHNE REICHWEITE C---- 20 N1 = K DX1 = EK - RGDATA(1,N1) RANGE = DX1**3 * RGDATA(5,N1) + 1 DX1**2 * RGDATA(4,N1) + 2 DX1 * RGDATA(3,N1) + 3 RGDATA(2,N1) RGRENZ = RANGE - RT IF ( RGRENZ .LE. 0.0D00 ) GOTO 15 C----- C SUCHE ENDINTERVALL C----- DO 30 K=N1,1,-1 IF ( RGDATA(2,K) .LE. RGRENZ ) GOTO 40 30 CONTINUE c-kt die Fehlerbehandlung wird der aufrufenden Routine überlassen: WRITE(SYSOUT,3007) 3007 FORMAT(//' Schwerer Fehler in der Funktion ''ELOSS2'':') WRITE(SYSOUT,3008) 3008 FORMAT(' ',T5,'Ein Teilchen hatte einen kleineren Wert für die'// 1 'verbleibende Reichweite'/T5,'als der kleinste tabellierte '// 2 'Reichweiten-Wert.') WRITE(SYSOUT,3009) RGRENZ, RGDATA(2,K) 3009 FORMAT(' ',T5,'verbleibende Reichweite des Teilchens',T50, 1D22.16,/T5,'kleinster tabellierter Reichweiten-Wert',T50,D22.16) ELOSS2 = -2. RETURN c-kt Ende c-kt STOP 11 C----- C ENDINTERVALL GEFUNDEN, SUCHE ENDENERGIE C----- 40 N1 = K C C NORMALFORM: X**3 + R*X**2 + S*X + T = 0 C R = RGDATA(4,N1) / RGDATA(5,N1) S = RGDATA(3,N1) / RGDATA(5,N1) T = ( RGDATA(2,N1) - RGRENZ ) / RGDATA(5,N1) C C REDUZIERTE DARSTELLUNG: Y**3 + P*Y + Q = 0; X = Y -(R/3) C c-kt P = ( 3.0D00 * S - R*R ) * DRITTL c-kt ersetzt durch P = S - R*R * DRITTL PDRITT = P * DRITTL RDRITT = R * DRITTL Q = (2.0D00 * RDRITT**3) - (RDRITT*S) + T D = PDRITT**3 + (0.25D00*Q**2) IF ( D .LT. 0.0D00 ) GOTO 50 C C EINE REELLE LOESUNG C U = ( -0.5D00*Q + DSQRT(D) ) ** DRITTL V = -PDRITT / U EKOUT = U + V - RDRITT + RGDATA(1,N1) GOTO 70 C C DREI REELLE LOESUNGEN C 50 RHO = DSQRT ( -PDRITT**3 ) PHI = DACOS ( -Q*0.5D00/RHO ) RHO = 2.0D00 * RHO**DRITTL C XDUM = RGDATA(1,N1+1) - RGDATA(1,N1) c-kt Die Genauigkeit von 4*PI/3 ist von 8 auf 16 Stellen erweitert: X = RHO * DCOS ( (PHI/3.0D00) + 4.188790204786391D00 ) - RDRITT IF ( X .GE. 0.0D00 .AND. X .LE. XDUM ) GOTO 60 X = RHO * DCOS ( PHI/3.0D00 ) - RDRITT IF ( X .GE. 0.0D00 .AND. X .LE. XDUM ) GOTO 60 c-kt Die Genauigkeit von 2*PI/3 ist von 8 auf 16 Stellen erweitert: X = RHO * DCOS ( (PHI/3.0D00) + 2.094395102393195D00 ) - RDRITT IF ( X .GE. 0.0D00 .AND. X .LE. XDUM ) GOTO 60 c-kt die Fehlerbehandlung wird der aufrufenden Routine überlassen: WRITE(SYSOUT,3010) 3010 FORMAT(/' Schwerer Fehler in der Funktion ''ELOSS2'':') WRITE(SYSOUT,3011) 3011 FORMAT(' ',T5,'Durch Rundungsfehler wurde keine korrekte Lösung' 1/T5,'des kubischen Spline-Polynoms gefunden.') ELOSS2 = -3. RETURN c-kt Ende c-kt STOP 9 C 60 EKOUT = X + RGDATA(1,N1) C 70 ELOSS2 = EK - EKOUT RETURN END C%CREATE INDEDX NOKEYS C----------------------------------------------------------------------- C C NAME: INDEDX C C FUNKTION: BERECHNET 1/DEDX-FUNKTION FUER BELIEBIGE GEMISCHE C C AUFRUF: VDEDX = INDEDX ( EKIN ) C C PARAMETER: C C EKIN R*8 KIN. ENERGIE DES TEILCHENS IN MEV C C FUNKTIONSWERT: C C VDEDX R*4 1. / DEDX C C BEMERKUNG: C C (1) ES WIRD DIE FUNKTION "DEDXM" AUFGERUFEN, D.H. ES MUESSEN C VORHER ALLE PARAMETER AUF DEM COMMON-BLOCK /DEDXPA/ C INITIALISIERT WERDEN. C C (2) DIE DATEN MUESSEN AUF DEM COMMON-BLOCK /DEDXPA/ C WIE FOLGT ABGELEGT WORDEN SEIN: C C SMASS R*4 MASSE DER TEILCHENS IN MEV/C**2 C NZ I*4 KERNLADUNGSZAHL (INTERN FUER DEDXM!) C NEL I*4 ZAHL DER ELEMENTE C WMOL R*4 "MOLEKULARGEWICHT" = SUMME ( ATOM(I)*AWT(I) ) C NZL R*4 KERNLADUNGSZAHLEN DER ELEMENTE C ATOM R*4 ZAHL DER ATOME DER ELEMENTE C C----------------------------------------------------------------------- REAL*8 FUNCTION INDEDX(EKIN) real*8 ekin EXTERNAL DEDXM C---- INDEDX = 0.0d+00 DUMMY = DEDXM(sngl(EKIN)) IF ( DUMMY .EQ. 0.0 ) GOTO 10 INDEDX = 1.0 / DUMMY 10 RETURN END C%CREATE MLR SUBROUTINE MLR(Z,AV,RHA,TM,EN,PMASS,ZE,FI) EXTERNAL GGUBS INTEGER*4 SEEDS, TEST c-kt Deklaration eines neuen Bezeichners für die Ausgabeunit: SYSOUT INTEGER*2 SYSOUT, SYSIN COMMON 1 / STATUS / SYSOUT, SYSIN, TEST COMMON / SEEDCO / SEEDS(6) C--------------------------------------------------------------------- C C ALLE DATA-KONSTANTEN MIT E-EXPONENTEN VERSEHEN. C S.VETTER URZ HEIDELBERG 8.2.1980 C C GENERATION OF MULTIPLE SCATTERING (MOLIERE FORMULA) C *************************************************** C Z= ATOMIC NUMBER OF SCATTERER C AV= ATOMIC WEIGHT OF THE SCATTERER C RHA= DENSITY OF THE SCATTERER C TM= SCATTERER THICKNESS (CM) C EN= INCIDENT PARTICLE MOMENTUM (MEV/C) C PMASS= MASS OF INCIDENT PARTICLE (MEV) C ZE= CHARGE OF INCIDENT PARTICLE (IN UNITS OF THE ELECTRON CHARGE) C FI= GENERATED SCATTERING ANGLE C ATA(NA),FK1(NA),FK2(NA),FK3(NA)= TABLES USED IN MOLIERE FORMULA C FOR MULTIPLE SCATTERING C C--------------------------------------------------------------------- DIMENSION A(60),ARG(60),VAL(4),ATA(55),FK1(55),FK2(55),FK3(55) DATA ATA/0.0, .0999999 , .1999999 , .2999998 , .3999999 , * .4999997 , .5999999 , .6999997 , .7999998 , .8999996 , .9999998 , *1.199999 ,1.399999 ,1.599999 ,1.799999 ,1.999999 ,2.199999 , *2.399999 ,2.599999 ,2.799999 ,2.999999 ,3.199999 ,3.399999 , *3.599999 ,3.799999 ,3.999999 ,4.199998 ,4.399999 ,4.599999 , *4.799999 ,4.999999 ,5.199998 ,5.399999 ,5.599999 ,5.799999 , *5.999999 ,6.199998 ,6.399999 ,6.599999 ,6.799999 ,6.999999 , *7.199998 ,7.399997 ,7.599999 ,7.799999 ,7.999999 ,8.199998 , *8.399997 ,8.599999 ,8.799999 ,8.999999 ,9.199998 ,9.399997 , *9.599999 ,9.799999 / DATA FK1/ .0, .00995016E0, .03921056E0, .08606875E0, .1478561 , * .2211991 , .3023236 , .3873733 , .4727075 , .5551416 , .6321204 , * .7630718 , .8591412 , .9226950 , .9608359 , .9816843 , .9920929 , * .9968488 , .9988408 , .9996063 , .9998766 , .9999643 , .9999904 , * .9999976 , .9999994 , .9999998 ,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1., *1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1./ DATA FK2/0.0, .00415007E0, .01550778E0, .03102755E0, .04642320E0, * .05693358E0, .05814465E0, .04697935E0, .02206861E0,- .01593698E0, *- .0645485 ,- .1776276 ,- .2818608 ,- .3497007 ,- .3716551 , *- .3548789 ,- .3147738 ,- .2664943 ,- .2203767 ,- .1812971 , *- .1502245 ,- .1261782 ,- .1076097 ,- .0930889 ,- .0815174 , *- .0721195 ,- .0643684 ,- .0578844 ,- .0523495 ,- .0476232 , *- .0435265 ,- .0399468 ,- .0368066 ,- .0340244 ,- .0315585 , *- .0293560 ,- .0273818 ,- .0256116 ,- .0240157 ,- .0225666 , *- .0212406 ,- .0200279 ,- .0189221 ,- .0179088 ,- .0169753 , *- .0161105 ,- .0153094 ,- .0145689 ,- .0138822 ,- .0132431 , *- .0126462 ,- .0120869 ,- .0115616 ,- .0110680 ,- .0106048 / DATA FK3/ .0E0, .01223041E0, .04564741E0, .09133768E0, .1373360E0, *.1715474 ,.1842774 ,.1711607 ,.1335052 ,.07785457E0,.01460415E0, *-.09175974E0,-.1223490E0,-.07293552E0,.01300466E0,.08469558E0, *.1153594,.1071076,.07755226E0,.04448753E0,.01822135E0,.00146502E0, *-.734811E-2,-.0109819,-.1178313E-1,-.1125726E-1,-.1019054E-1, *-.896781E-2,-.786496E-2,-.683298E-2,-.592750E-2,-.514682E-2, *-.446538E-2,-.387442E-2,-.336000E-2,-.291120E-2,-.251946E-2, *-.217772E-2,-.187836E-2,-.161443E-2,-.137980E-2,-.117159E-2, *-.098778E-2,-.082460E-2,-.067865E-2,-.054704E-2,-.042855E-2, *-.032243E-2,-.022696E-2,-.014061E-2,-.006204E-2,.000983E-2, *.007585E-2,.013652E-2,.019206E-2/ BETA=EN/SQRT(PMASS**2+EN**2) IF(PMASS.LT.1.0) Q=SQRT(Z*(Z+1.)) IF(PMASS.GE.1.0) Q=Z RTN=SQRT(RHA*TM/AV) IF(PMASS.LT.1.0) QSML=(Z+1.)*Z**(1./3.) IF(PMASS.GE.1.0) QSML=Z**(4./3.) c ( Molière (CHIa/CHI0)**2 ) = 'EPSN' c ( Molière (SHIc/CHIa)**2 ) = 'GAMMA' c ( Molière b ) = 'CNST' c EPSN=1.13+3.76*(Z*ZE/(137.*BETA))**2 c GAMMA=(8.831E3*QSML*RHA*TM)/(BETA**2*AV*EPSN) c CNST=ALOG(GAMMA)-0.154 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc EPSN = 1.13 + 3.76 * ( Z*ZE / (137.036*BETA) )**2 GAMMA = 8838.3647*QSML*RHA*TM/(BETA**2*AV*EPSN) CNST = ALOG(GAMMA) - 0.15443133 IF ( CNST .LE. 1.0 ) THEN WRITE(SYSOUT,1) TM, RHA, Z, CNST, 1.167*EXP(CNST), 3.2/GAMMA*TM 1 FORMAT(' ',/' ',/' ***** UNZULÄSSIGER WERT IN DER '// 1 'EINGABEDATEI *****',/' ',/' Die Schicht mit der Dicke ', 2 E9.3,' cm, der Dichte ',E9.3,' g/cm**3 und der',/' mittleren '// 3 'Ladungszahl ',E9.3,' ist zu dünn, als daß die Näherungen '// 4 'des Modells,',/' das diesem Programm zugrundeliegt, noch '// 5 'Gültigkeit haben.',/' ',/' Die Molière-Konstante b beträgt ' 6 ,F8.7,'; die mittlere Anzahl',/' von Einzelstreuprozessen in '// 7 'der Schicht beträgt ',F6.4,'.',/' ',/' Das Modell wird '// 8 'rechenbar erst ab mindestens 3.2 Einzelstreuprozessen pro', 9 /' Schicht, dies entspricht einer Konstanten b von mindestens '// 1 '1.0.',/' Für die gegebenen Werte für Ladungszahl und Dichte '// 2 'müßte die Schichtdicke',/' auf ',E9.3,' cm erhöht werden, um '// 3 'sie in den berechenbaren Bereich zu bringen.',/' ',/ 4 ' ***** PROGRAMM ABGEBROCHEN *****') STOP 12 ENDIF c Lösung der transzendenten Gleichung B-ln(B)=b in max. 10 Schritten c (Molière B) = BM BM=5. L=0 10 DBM=-(BM-ALOG(ABS(BM))-CNST)/(1.-1./BM) BM=BM+DBM L=L+1 IF(L-10) 20,20,190 20 IF(ABS(DBM)-0.0001) 30,30,10 30 FI1=3.965E-1*Q*ZE/(EN*BETA)*RTN CALL GGUBS ( SEEDS(1),1,YFL) IF(YFL.GT.(1.+FK2(55)/BM)) GO TO 180 A(4)=FK1(4)+FK2(4)/BM+FK3(4)/BM**2 NA=0 40 NA=NA+1 A(NA)=FK1(NA)+FK2(NA)/BM+FK3(NA)/BM**2 S=YFL-A(NA) IF(S) 60,60,50 50 CONTINUE GO TO 40 60 CONTINUE IF(NA.EQ.55) NA=54 A(NA+1)=FK1(NA+1)+FK2(NA+1)/BM+FK3(NA+1)/BM**2 C C INTERPOLATION ACCORDING TO AITKEN"S METHOD C ****************************************** C YFL= ARGUMENT VALUE C ARG= ARGUMENT VALUES OF THE TABLE (ARG,VAL) C VAL= FUNCTION VALUES OF THE TABLE (ARG,VAL) C NDIM= NUMBER OF POINTS IN TABLE (ARG,VAL) C ATI= THE RESULTING INTERPOLATED FUNCTION VALUE C NDIM=4 IF(NA.LT.3) NA=3 DO 70 M=1,NDIM ARG(M)=A(NA-NDIM+1+M) VAL(M)=ATA(NA-NDIM+1+M) 70 CONTINUE F=ATA(NA)/100. DELT2=0. 80 DO 130 J=2,NDIM DELT1=DELT2 IEND=J-1 DO 90 I=1,IEND H=ARG(I)-ARG(J) IF(H) 90,140,90 90 VAL(J)=(VAL(I)*(YFL-ARG(J))-VAL(J)*(YFL-ARG(I)))/H DELT2=ABS(VAL(J)-VAL(IEND)) IF(J-2) 130,130,100 100 IF(DELT2-F) 160,160,110 110 IF(J-5) 130,120,120 120 IF(DELT2-DELT1) 130,140,140 130 CONTINUE GO TO 150 140 J=IEND GO TO 160 150 J=NDIM 160 ATI=VAL(J) 170 FI=ATI*FI1*SQRT(BM) GO TO 200 180 ATI=SQRT(5./(1.-(1.-BM*(1.-YFL))**5)) FI=ATI*FI1*SQRT(BM) GO TO 200 190 FI=0. 200 RETURN END C%CREATE STRAG1 C----------------------------------------------------------------------- C C NAME: STRAG1 C C FUNKTION: BERECHNET DAS RELATIVE RANGE-STRAGGLING C C PARAMETER: P R*4 IMPULS IN MEV/C C EM R*4 MASSE IN MEV/C**2 C C BEMERKUNG: DIESES PROGRAMM IST EINE VERAENDERTE VERSION DES C CERNLIB-PROGRAMMS "STRAG". DER NAME WURDE AUS C SICHERHEITSGRRUENDEN IN "STRAG1" GEAENDERT, DIE C PARAMETERLISTE WURDE MEHR DEN ERFORDERNISSEN DER C "MITTELENERGIEPHYSIK" ANGEPASST. C C----------------------------------------------------------------------- FUNCTION STRAG1(P,EM) DIMENSION ENT(15),TAB(15) DATA ENT/.0,.001,.00189,.00355,.00675,.0128,.0282,.0460,.0875,.16 + 5,.312,.593,1.012,5.,100000./ DATA TAB/.1,.01573,.01445,.0134,.0124,.0117,.0112,.0107,.0102,.00 + 97,.0092,.0087,.0082,0.,0./ BI = SQRT(1.+(P/EM)**2)-1. DO 1 I=2,15 IF ( ENT(I) .GT. BI ) GO TO 2 1 CONTINUE 2 A = (BI-ENT(I-1))*(TAB(I)-TAB(I-1))/(ENT(I)-ENT(I-1))+TAB(I-1) STRAG1 = SQRT( 938.213/EM)*A RETURN END C----------------------------------------------------------- C C RANDOM NUMBER GENERATOR FOR COMPATIBILITY C C WITH IBM ROUTINES C C UNIFORM DISTRIBUTION C C CALL GGUBS(K,N,R) C C IN R(N) ARE THE RANDOM NUMBERS C C----------------------------------------------------------- SUBROUTINE GGUBS(K,N,R) REAL R(N) c DO 100 I=1,N c R(I)=RAN(K) c100 CONTINUE call ranlux(R,N) RETURN END C--------------------------------------------------------- C C GAUSSIAN DISTRIBUTION C C CALL GGNML (K,N,R) C C N <=2 C C IN R(N) ARE THE NUMBERS STORED C--------------------------------------------------------- SUBROUTINE GGNML(K,N,R) REAL R(N) real random(2) 10 call ranlux(random,2) A = random(1) c10 A=RAN(K) IF(A.LE.0.) GO TO 10 c B=RAN(K) B = random(2) C=SQRT(-2.*ALOG(A)) D=6.2831853*B R(1)=C*COS(D) IF(N.EQ.2) R(2)=C*SIN(D) RETURN END