Files
musrsim/mutrack/src/SUB_ELOSS.FOR
2005-03-22 10:33:08 +00:00

664 lines
19 KiB
Fortran

c===============================================================================
OPTIONS /EXTEND_SOURCE
SUBROUTINE CALC_ELOSS_ICRU(E0,q,m,Thickness,Everlust)
c =====================================================
c Einheiten:
c ----------
c Energie: [E0] = keV
c Ladung: [q] = e
c Masse: [m] = keV/c**2, wird auf m(Proton) umgerechnet
c Foliendicke [Thickness] = ug/cm**2 ! => 'masse'
c Energieverlust [Everlust] = keV
IMPLICIT NONE
real E0,EVerlust ! Startenergie und Energieverlust
real E,q,m,masse,Thickness ! Energie,Masse,Ladung,Foliendicke
real x ! aktuelle Position in der Folie
real dx ! Vorgabe fuer Dickenschritt
COMMON /ELoss_Ex/ E,x
COMMON /ELoss_masse / masse
integer steps ! Nummer des aktuellen Integrationsschrittes
integer maxSteps /2000/ ! maximal tolerierte Anzahl an Schritten
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c Ueberpruefen ob (q.EQ.+1), sonst macht diese Energieverlustberechnung keinen
c Sinn:
if (q.NE.1) then
write(*,*)
write(*,*) 'SUBROUTINE CALC_ELOSS_ICRU: Teilchenladung q = ',q
write(*,*)
write(*,*) 'Berechnung des Energieverlustes aus Protonendaten nur'
write(*,*) 'fuer Teilchen mit Ladungszahl q = +1 implementiert'
write(*,*)
write(*,*) '-> STOP'
write(*,*)
STOP
endif
c Startwerte setzen:
masse = m / 938272.3 ! Umrechnen in Protonenmassen
E = E0 ! aktuelle Energie auf Startenergie setzen
x = 0. ! Wir beginnen auf der Folienforderseite
dx = .05 ! Dring beim ersten Versuch 0.05 Mikrogramm/cm**2 in
! die Folie ein
steps = 0. ! Zaehler fuer Anzahl der Integrationsschritte resetten
c...............................................................................
c hierher wird zurueckgesprungen, solange Folienende noch nicht erreicht ist.
10 if (x+dx.GT.Thickness) dx = Thickness-x ! Schritt erreicht Folienende
Steps = Steps + 1
call INTEGRATION_ELOSS(dx)
if (x.EQ.Thickness) then
EVerlust = E0 - E
c write(*,*) 'Ein, ELoss = ',E0,EVerlust
c write(*,'(A,i4,A)') ' calculation of ''ELOSS'' needed ',steps,' steps'
RETURN
elseif (x.GT.Thickness+1.E-5) then
EVerlust = E0 - E
write(*,*)
write(*,*) 'SUBROUTINE CALC_ELOSS_ICRU: x > Foliendicke'
write(*,*)
write(*,*) 'x = ',x
write(*,*) 'Thickness = ',Thickness
write(*,*) 'x-Thickness = ',x-Thickness
write(*,*) 'E0 = ',E0
write(*,*) 'E = ',E
write(*,*) 'EVerlust = ',EVerlust
write(*,*) 'maxSteps = ',maxSteps
write(*,*) 'steps = ',steps
write(*,*)
write(*,*) '-> STOP'
write(*,*)
STOP
elseif (Steps.GT.MaxSteps) then
write(*,*)
write(*,*) 'SUBROUTINE CALC_ELOSS_ICRU:'
write(*,*)
write(*,*) 'Berechnung des Energieverlustes nicht innerhalb der erlaubten'
write(*,*) 'Maximalzahl an Schritten vollendet:'
write(*,*)
write(*,*) 'steps = ',steps
write(*,*) 'maxSteps = ',maxSteps
write(*,*) 'x = ',x
write(*,*) 'Thickness = ',Thickness
write(*,*) 'Thickness-x = ',Thickness-x
write(*,*)
write(*,*) '-> STOP'
write(*,*)
STOP
else
goto 10
endif
END
c===============================================================================
OPTIONS /EXTEND_SOURCE
SUBROUTINE INTEGRATION_ELOSS(dx)
c ================================
IMPLICIT NONE
SAVE
c Diese Subroutine berechnet zu einem vorgegebenen Dickenschritt dx den
c Energieverlust zweimal: einmal direkt mit dx und einmal ueber zwei
c aufeinanderfolgende Schritte mit dx/2. (die beiden dx/2-Schritte werden
c zuerst ausgefuehrt).
c
c Aus der Differenz der beiden Resultate wird eine Abschaetzung fuer den Fehler
c des dx-Schrittes gewonnen, die dazu verwendet wird zu entscheiden, ob der
c Integrationsschritt mit einem verkuerzten Dickenintervall wiederholt werden
c muss, oder ob das Dickenintervall fuer den folgenden ausgedehnt werden kann.
c
c Die beiden Einzelergebnisse aus dem dx- und den beiden dx/2-Schritten, die
c jeweils ueber Runge-Kutta-Rechnung vierter Ordnung erhalten werden, werden
c zum Schluss noch zusammengenommen, um ein Resultat mit Genauigkeit fuenfter
c Ordnung in dx zu erhalten.
c
c Der ganze Ablauf entspricht den Ausfuehrungen in Kapitel 15.2 der NUMERICAL
c RECIPIES: 'Adaptive Stepsize Control for Runge-Kutta'.
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
real E ! aktuelle Energie bei Aufruf und
! Verlassen der Routine
real x ! aktuelle Position in der Folie
COMMON /ELoss_Ex/ E,x
logical log_relativ /.false./ ! relative Fehlerbetrachtung?
real eps /1e-8/ ! verlangte Genauigkeit des einzelnen
! Integrationsschrittes
real dxsmall /.001/ ! minimale Dickenaenderung
integer n_dxSmall ! wie oft wurde dxSmall unterschritten
integer maxBelowDxSmall / 50 / ! max. tolerierte Anzahl an Unterschr.
real dx ! Dickenschritt
real dx_half ! halber Dickenschritt
real dE_dx_0,dE_dx_1 ! Stoppingpower
real E1,DE1,DE2 ! fuer Energieintegration
real EDifferenz ! fuer Fehlerbetrachtung
real Error ! fuer Fehlerbetrachtung
real errCon, safety ! fuer Schrittweitensteuerung
real pShrink, pGrow ! fuer Schrittweitensteuerung
PARAMETER (errCon = 6.e-4, safety = .9) ! vgl. Referenz
PARAMETER (pShrink = -.25, pGrow = -.2)
! errCon = (4./safety)**(1/pGrow)
logical flag_dxSmall ! wenn bei dx < dxsmall der Fehler
! immer noch zu gross ist.
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
flag_dxSmall = .false. ! flag resetten
c Lese stopping power bei aktueller Energie. Speichere sie in dE_dx_0, damit sie
c wiederverwendet werden kann, falls die Berechnung mit kuerzerem Dickenschritt
c wiederholt werden muss.
call stoppingPower(E,dE_dx_0)
c...............................................................................
10 continue ! hier gehts wieder von vorne los, falls Dickenschritt dx
! abgeaendert werden muss.
dx_half = dx / 2.
c mache ersten dx/2 - Schritt:
call RK4_ELOSS(dx_half,E,dE_dx_0,DE1)
c Lese stopping power bei Energie E1=E+DE1 und mache zweiten dx/2 - Schritt:
E1 = E+DE1
call stoppingPower(E1,dE_dx_1)
call RK4_ELOSS(dx_half,E1,dE_dx_1,DE2)
c Summiere Energieverluste der beiden dx/2 -Schritte und speichere in DE1:
DE1 = DE1 + DE2
c mache dx - Schritt:
call RK4_ELOSS(dx,E,dE_dx_0,DE2)
c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c Fehlerbetrachtung und Berechnung des endgueltigen Ergebnisses:
C Fehlerbetrachtung: der groesste (absolute bzw. relative) Fehler soll kleiner
C als eps sein:
Error = 0.
EDifferenz = DE1-DE2
if (log_relativ) then
if (DE1.NE.0.) then
Error = Abs( EDifferenz/DE1 )
elseif (DE2.NE.0.) then
Error = Abs( EDifferenz/DE2 )
endif
else
Error = Abs( EDifferenz )
endif
c - Skaliere Fehler auf Epsilon:
Error = Error / eps
c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c Die Groesse des reskalierten 'Error' bestimmt, ob der Integrationsschritt
c mit kleinerem Dickenintervall wiederholt werden muss, bzw. um welchen
c Faktor das Dickenintervall fuer den naechsten Schritt vergroessert werden
c kann:
c Liegt der Fehler ausserhalb des Toleranzbereiches und ist dx bereits jetzt
c kleiner als dxsmall, so mache keinen neuen Versuch sondern akzeptiere als Not-
c loesung den bestehenden Naeherungswert. Setze dx in diesem Fall als Default
c fuer den kommenden Integrationsschritt auf dxsmall. Setze aber auch das flag
c 'flag_dxsmall', damit gezaehlt werden kann, wie oft dieses Prozedur angewendet
c werden muss. Ist dies zu oft der Fall, so brich die Berechnung ganz ab.
c (2. Teil erfolgt weiter unten)
if (Error.GT.1.) then
if (dx.LT.dxsmall) then ! Fehler immer noch zu gross, obwohl
flag_dxsmall = .true. ! dxsmall schon unterschritten ist
else
!c Bestimme kuerzeren Dickenschritt fuer neuen Versuch (vgl. Referenz):
dx = safety * dx * (Error**pShrink)
goto 10
endif
endif
c Nimm die Ergebnisse aus dem dx-Schritt und den beiden dx/2-Schritten und
c berechne damit den Energieverlust mit Genauigkeit fuenfter Ordnung in dx:
E = E + DE1 + EDifferenz / 15.
c alten Dickenschritt addieren, neuen so gross wie sinnvoller weise moeglich
c machen:
x = x + dx
if (flag_dxSmall) then
if (n_dxsmall.LT.maxBelowdxSmall) then
dx = dxSmall ! fuer naechsten RK-Schritt
n_dxsmall = n_dxsmall + 1
else
write(*,*) 'n_dxSmall ueberschritten! -> STOP'
STOP
endif
else
if (Error.GT.errCon) then
dx = safety * dx * (Error**pGrow) ! vgl. Referenz
else
dx = 4. * dx ! <- Vergroesserung des Dickenschritts max. um
endif ! Faktor 4!
endif
END
c===============================================================================
OPTIONS /EXTEND_SOURCE
SUBROUTINE RK4_ELOSS(dx_,EIn,dE_dx_0,DE)
c ========================================
IMPLICIT NONE
SAVE
c Diese Subroutine berechnet bei vorgegebenem Dickenschritt einen einzelnen
c Runge-Kutta-Integrationsschritt 4. Ordnung.
c Die Vorgehensweise entspricht den Ausfuehrungen in Kapitel 15.1 der
c NUMERICAL RECIPIES: 'Runge-Kutta Method'. Die Formeln sind zwar teilweise
c etwas umgeschrieben, sind aber mathematisch alle aequivalent zu denen der
c Referenz.
c Zurueckgegeben wird die errechnete Energieaenderung anstatt direkt des neuen
c Wertes, da sonst Schwierigkeiten auftreten koennen, wenn in der Subroutine
c 'INTEGRATIONSSTEP_ELOSS' aus der Differenz der neuen Werte aus den beiden
c dx/2- und dem dx-Schritt der Fehler abgeschaetzt werden soll (kleine
c Differenz moeglicherweise grosser Werte).
real dE_dx_0, EIn ! Eingangsgroessen
real dx_,dx_half_ ! Dickenschritt, halber Dickenschritt
real DE ! Ergebnisspeicher
real Test_E ! Test-Energie
real dE_dx_1,dE_dx_2,dE_dx_3 ! Stoppingpower bei Test-Energien
c First step of Runge-Kutta-Method:
dx_half_ = dx_ / 2.
Test_E = EIn + dE_dx_0 * dx_half_
c Second step of Runge-Kutta-Method:
call stoppingPower(Test_E,dE_dx_1)
Test_E = EIn + dE_dx_1 * dx_half_
c Third step of Runge-Kutta-Method:
call stoppingPower(Test_E,dE_dx_2)
Test_E = EIn + dE_dx_2 * dx_
c Fourth step of Runge-Kutta-Method:
call stoppingPower(Test_E,dE_dx_3)
c calculate forward step:
DE = (dE_dx_0 + 2.*(dE_dx_1+dE_dx_2) + dE_dx_3) * dx_/6.
END
c===============================================================================
OPTIONS /EXTEND_SOURCE
SUBROUTINE stoppingPower(E_,dE_dx)
c =================================
IMPLICIT NONE
real E_,dE_dx,StoPo_
integer DatenZahl
parameter (DatenZahl=56)
real E(DatenZahl)
real StoPoGr(DatenZahl) ! fuer Graphit
real StoPoAK(DatenZahl) ! fuer amorphen Kohlenstoff
real Ep ! Energie von Protonen mit gleicher Geschwindigkeit
real masse ! Masse der Teilchen in Protonenmassen
COMMON /ELoss_masse / masse
integer i
logical graphitData
common /elossDataType/ graphitData
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c Zu den Einzelbeitraegen von Elektronen und Kernen vergleiche die Programme
c NUCL_AND_EL_STOP_AMORPH.FOR und NUCL_AND_EL_STOP_GRAPHIT.FOR im directory
c UD0:[SIMULA.PROGRAMS.ICRU-ELOSS] auf der SLOMU in Konstanz.
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
c Die Daten sind aus dem ICRU Report 49 entnommen:
c
c ICRU Report 49: 'Stopping Powers and Ranges for Protons and Alpha Particles'
c Signatur PSI-Bibliothek West: 850 B 1:49
c Tabelle Seite 112: PROTONS IN CARBON (Graphite, density 1.7 g/cm^3)
c Spalte: TOTAL STOPPING POWER (= electronic plus nuclear)
DATA E / 0.,
+ 001,001.5,002,002.5,003,004,005,006,007,008,009,
+ 010,012.5,015,017.5,020,022.5,025,027.5,030,035,040,
+ 045,050,055,060,065,070,075,080,085,090,095,
+ 100,125,150,175,200,225,250,275,300,350,400,
+ 450,500,550,600,650,700,750,800,850,900,950 /
c Die hier aufgelisteten Werte sind um Faktor 10000. zu gross. -> Wird unten
c reskaliert.
c Dies sind die totalen Stoppingpowerdaten fuer Graphit!
c Diese wurden bis Version 1.5.6 von Mutrack ausschliesslich verwendet.
DATA StoPoGr / 0.,
+ 1677,1963,2214,2443,2652,2958,3241,3500,3738,3946,4132,
+ 4303,4678,4991,5264,5510,5731,5927,6101,6259,6551,6823,
+ 7067,7274,7441,7574,7677,7755,7810,7845,7861,7860,7846,
+ 7819,7548,7153,6721,6297,5908,5561,5254,4982,4543,4196,
+ 3906,3654,3433,3240,3073,2925,2793,2675,2568,2470,2380 /
c Dies sind die totalen Stoppingpowerdaten fuer amorphen Kohlenstoff!
DATA StoPoAK / 0.,
+ 1635,1906,2140,2351,2544,2826,3086,3323,3541,3730,3899,
+ 4055,4393,4675,4919,5140,5338,5512,5666,5806,6066,6307,
+ 6523,6706,6854,6970,7061,7129,7178,7208,7221,7222,7209,
+ 7185,6955,6619,6254,5897,5566,5265,4996,4755,4362,4048,
+ 3782,3549,3343,3162,3004,2864,2739,2626,2524,2429,2343 /
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if (E_.LT.0) then
write(*,*)
write(*,*) 'error in subroutine ''STOPPINGPOWER'':'
write(*,*) 'E_ = ',E_,' < 0!'
write(*,*) '-> STOP'
write(*,*)
STOP
endif
c Energie auf diejenige Energie reskalieren, die ein Proton haette, wenn es
c sich mit der gleichen Geschwindigkeit bewegen wuerde, wie das Teilchen der
c Masse 'masse', das sich mit der Energie 'E_' bewegt:
Ep = E_/masse ! masse ist in Protonenmassen angegeben
c suche zustaendige Stelle in der Tabelle:
do i = 1, DatenZahl
if (E(i).GE.Ep) goto 10
enddo
write(*,*)
write(*,*) 'error in subroutine ''STOPPINGPOWER'': Ep > Ep_max'
write(*,*) '-> STOP'
write(*,*)
STOP
10 continue
c Berechne die Stoppingpower (ueber lineare Interpolation):
if (i.lt.2) then
write(*,*) 'i = ',i
STOP
endif
if (graphitData) then
StoPo_ = StoPoGr(i-1) + (StoPoGr(i)-StoPoGr(i-1))*(Ep-E(i-1))/(E(i)-E(i-1))
else
StoPo_ = StoPoAK(i-1) + (StoPoAK(i)-StoPoAK(i-1))*(Ep-E(i-1))/(E(i)-E(i-1))
endif
dE_dx = - StoPo_ / 10000. ! negativieren und reskalieren
END
c===============================================================================
OPTIONS /EXTEND_SOURCE
SUBROUTINE E_Straggling_Yang(mean_E,m,sigmaEloss)
c =================================================
IMPLICIT NONE
c Dieses Routine berechnet die Energieverluststreuung fuer einfach geladenen
c Teilchen in der Kohlenstofffolie entsprechend den empirischen Formeln von
c Q. Yang: "Empirical formulae for energy loss straggling of ions in matter",
c Q. Yang, D.J. O'Connor, Y. Wang, NIM B61, S. 149-155 (1991)
c Eingabeparameter: Foliendicke, Teilchenart und Eingangsenergie.
c Ausgabegroesse: Energieverluststraggling
c Einheiten:
c ----------
c Energie: [E0] = keV
c Ladung: [q] = e
c Masse: [m] = keV/c**2, wird auf m(Proton) umgerechnet
c Foliendicke [Thickness] = ug/cm**2 ! => 'masse'
c Straggling: [sigmaEloss] = keV
INCLUDE 'mutrack$sourcedirectory:COM_MUTRACK.INC'
real mean_E ! mittlere Energie in Folie
real m ! Masse
real mean_Ep ! mean_E umgerechnet auf Protonenmasse
! (fuer sigmaE = 0)
real sigmaEloss ! Energieverluststreuung (Ausgabevariable)
c real sigmaInhom ! Streuung aufgrund von Dickeninhomog.
real Z2
parameter (Z2 = 6.) ! Kernladungszahl von Kohlenstoff
real factor_Thickness ! fuer Umrechnung von ug/cm**2
parameter (factor_Thickness=6.022e23/12.011e6) ! in Atome/cm**2
! 6.022e23: Avogadrokonstante
! 12.011e6: molare Masse in ug.
real eSquare ! fuer Umrechnung in SI-System
parameter (eSquare = 1.44e-10) ! keV*cm
real fourPi
parameter (fourPi = 12.56637061)
real factor_Bohr ! fuer Berechnung der Streuung nach Bohr
parameter (factor_Bohr = fourPi*Z2*eSquare*eSquare*factor_Thickness)
real OmQuad, OmQuadBohr
real OmQuad_over_OmQuadBohr
real DOmQuad_over_OmQuadBohr
real OmQuadChu_over_OmQuadBohr
real Gamma, rHelp
real A1,A2,A3,A4
parameter (A1 = -0.5127, A2 = -0.8595, A3 = 0.5626, A4 = -0.8721)
real B1,B2,B3,B4
parameter (B1 = 0.1955, B2 = 0.6941, B3 = 2.522, B4 = 1.040)
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c Umrechnen der Energie auf Protonenenergie in MeV:
mean_Ep = mean_E * 938272.6/m / 1000. ! mean_Ep in MeV !!!
c Berechnen der Streuung:
! Streuung nach Bohr:
OmQuadBohr = factor_Bohr * Thickness
! Streuung nach Chu:
OmQuadChu_over_OmQuadBohr = 1. / (1. + A1*mean_Ep**A2 + A3*mean_Ep**A4)
! Streuung nach Yang:
Gamma = B3 * (1. - exp(-B4*mean_Ep))
rhelp = mean_Ep - B2
DOmQuad_over_OmQuadBohr = B1*Gamma / (rHelp*rHelp + Gamma*Gamma)
! Gesamtresultat:
OmQuad_over_OmQuadBohr = OmQuadChu_over_OmQuadBohr + DOmQuad_over_OmQuadBohr
OmQuad = OmQuad_over_OmQuadBohr * OmQuadBohr
c Folieninhomogenitaeten beruecksichtigen:
c OmQuad = OmQuad + (sigmaInhom * sigmaInhom)
c Endgueltige Energiestreuung:
sigmaEloss = SQRT(OmQuad)
END
c===============================================================================
OPTIONS /EXTEND_SOURCE
SUBROUTINE E_Straggling_Lindhard(mean_E,m,sigmaEloss)
c =====================================================
IMPLICIT NONE
c Dieses Routine berechnet die Energieverluststreuung fuer einfach geladenen
c Teilchen in der Kohlenstofffolie entsprechend den Formeln von J. Lindhard
c und M. Scharff
c Eingabeparameter: Foliendicke, Teilchenart und Eingangsenergie.
c Ausgabegroessen: Energieverluststraggling
c Einheiten:
c ----------
c Energie: [E0] = keV
c Ladung: [q] = e
c Masse: [m] = keV/c**2, wird auf m(Proton) umgerechnet
c Foliendicke [Thickness] = ug/cm**2 ! => 'masse'
c Straggling: [sigmaEloss] = keV
INCLUDE 'mutrack$sourcedirectory:COM_MUTRACK.INC'
real mean_E ! mittlere Energie in Folie
real m ! Masse, Foliendicke
real sigmaEloss ! Energieverluststreuung (Ausgabevariable)
c real sigmaInhom ! Streuung aufgrund von Dickeninhomog.
real sigmaQuad
real help,L_y,y ! Fuer Berechnung der Lindhard-Scharff Korrektur
real v_over_v0_Quad
real fourPi
parameter (fourPi = 12.56637061)
real Z2
parameter (Z2 = 6.) ! Kernladungszahl von Kohlenstoff
real factor_v0Quad ! fuer Geschwindigkeitsskalierung auf v0
parameter (factor_v0Quad = 37574.632) ! (Bohrgeschwindigkeit)
! = 1/(13.6eV/me) , me in keV/c**2 eingesetzt
real factor_Thickness ! fuer Umrechnung von ug/cm**2
parameter (factor_Thickness=6.022e23/12.011e6) ! in Atome/cm**2
! 6.022e23: Avogadrokonstante
! 12.011e6: molare Masse in ug.
real eSquare ! fuer Umrechnung in SI-System
parameter (eSquare = 1.44e-10) ! keV*cm
real factor_Bohr ! fuer Berechnung der Streuung nach Bohr
parameter (factor_Bohr = fourPi*Z2*eSquare*eSquare*factor_Thickness)
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c Sigma nach Bohr berechnen:
sigmaQuad = factor_Bohr * Thickness
c Korrektur nach Lindhard & Scharff beruecksichtigen:
v_over_v0_Quad = mean_E/m * factor_v0Quad
y = v_over_v0_Quad / Z2
if (y.LT.3.) then
help = SQRT(y)
L_y = 1.36 * help - .016 *help*help*help
sigmaQuad = sigmaQuad * L_y / 2.
endif
c Folieninhomogenitaeten beruecksichtigen:
c sigmaQuad = sigmaQuad + (sigmaInhom * sigmaInhom)
c Endgueltige Energiestreuung berechnen:
sigmaEloss = SQRT(sigmaQuad)
END
c===============================================================================