unit unCross; interface const naMax=10000; type cross_ratio_type=array[1..naMax] of double; procedure get_cross_ratio(Z1_in,m1_in,Z2_in,m2_in,q_in,E_in,b_max_in : double; nb_in: integer; ep_in : double; theta_initial,dtheta : double; ntheta : integer; var cross_ratio : cross_ratio_type); implementation uses math; const nmax = 500000; type re = extended; var Z1,Z2 : re; ro : re; rmax : re; Ecm : re; b : double; ep : double; phi : array[0..nmax] of single; a_scr : re; Tmax, vmax : re; function V_pot(r : re) : re; var i : integer; begin if r > rmax then begin V_pot := 0; end else begin i := round(r/rmax*nmax); if i = nmax then V_pot := 0 else V_pot := Z1*Z2/r*exp((r-i*rmax/nmax)/(rmax/nmax) * (ln(phi[i+1])-ln(phi[i])) + ln(phi[i])); end; end; function V_pot_prime(r : re) : re; var i : integer; begin if r > rmax then begin V_pot_prime := 0; end else begin i := round(r/rmax*nmax); V_pot_prime := Z1*Z2*(-phi[i]/sqr(r) +(phi[i+1]-phi[i])/(r*rmax/nmax) ); end; end; procedure get_ro; var F,G : re; F_old : re; r, dr : re; count : integer; begin r := sqrt(sqr(b) + sqr(Z1*Z2/Ecm)); count := 0; F := 100; repeat inc(count); F_old := F; F := sqr(r)*(1 - V_pot(r)/Ecm) - sqr(b); G := 2*r*(1 - V_pot(r)/Ecm)-sqr(r)*V_pot_prime(r)/Ecm; dr := F/G; r := r - dr; if count > 1000 then begin writeln(count,' ',r,dr/r,F); end; until ((abs(dr/r) < 1e-14) and (count > 10))or ((count > 1000) and (F_old*F < 0)); // writeln(count); ro := r; end; function fun(x: re) : re; var r : re; begin r := ro/(1-sqr(x)); fun := x/sqrt(abs(1-V_pot(r)/Ecm-sqr(b/r))); end; function simp_inner(a1,b,ep : re): re; type arr = array[1..30] of re; var dx,epsp,x2,x3, f2,f3,f4,fmp, fbp,est2,est3 : arr; nrtr : array[1..30] of integer; pval : array[1..30,1..3] of re; sum,eps,fm,f1, a,absar, est,da, fa,fb,sx,est1 : re; lvl,l{,l1} : integer; label 1,2,3,4,5,7,11,12,13; begin { nonrecursive adaptive integration } { algorithm 182 CACM 6 (1983) 315 } { adaptation from : Numerical Integration , Davis-Rabinowitz, pg. 198 } a := a1; eps := ep; lvl := 0; absar := 0; est := 0; da := b-a; fa := 0;{fun(a);} fm := 4*fun((a+b)*0.5); fb := fun(b); { 1 = recur } 1: inc(lvl); dx[lvl] := da/3; sx := dx[lvl]/6; f1 := 4*fun(0.5*dx[lvl] + a ); x2[lvl] := a + dx[lvl]; f2[lvl] := fun(x2[lvl]); x3[lvl] := x2[lvl] + dx[lvl]; f3[lvl] := fun(x3[lvl]); epsp[lvl] := eps; f4[lvl] := 4*fun(dx[lvl]*0.5+x3[lvl]); fmp[lvl] := fm; est1 := sx*(fa+f1+f2[lvl]); fbp[lvl] := fb; est2[lvl] := sx*(f2[lvl] +f3[lvl] + fm); est3[lvl] := sx*(f3[lvl] + f4[lvl] + fb); sum := est1+ est2[lvl]+ est3[lvl]; absar :=absar - abs(est) + abs(est1) + abs(est2[lvl]) + abs(est3[lvl]); if abs(est-sum) - epsp[lvl]*absar > 0 then goto 3 else goto 2; 3: if lvl < 30 then goto 4 else goto 2; { 2=up } 2: dec(lvl); l := nrtr[lvl]; pval[lvl,l] := sum; case l of 1 : goto 11; 2 : goto 12; 3 : goto 13; else goto 4; end; 4: nrtr[lvl] := 1; est := est1; fm := f1; fb := f2[lvl]; 7: eps := epsp[lvl]/1.7; da := dx[lvl]; goto 1; 11: nrtr[lvl] := 2; fa := f2[lvl]; fm := fmp[lvl]; fb := f3[lvl]; est := est2[lvl]; a := x2[lvl]; goto 7; 12: nrtr[lvl] := 3; fa := f3[lvl]; fm := f4[lvl]; fb := fbp[lvl]; est := est3[lvl]; a := x3[lvl]; goto 7; 13: sum := pval[lvl,1] + pval[lvl,2] + pval[lvl,3]; if lvl>1 then goto 2 else goto 5; 5: simp_inner := sum; end; procedure get_theta(var theta : re); const r_cut : re = 4; var x_cut : re; begin r_cut := rmax; get_ro; r_cut := r_cut + ro; x_cut := sqrt(1-ro/r_cut); theta := pi - 4*b/ro*simp_inner(0,x_cut,ep/1000)-2*(pi/2-arccos(b/r_cut)); {theta := pi - 4*b/ro*romb(0,x_cut,ep/10)-2*(pi/2-arccos(b/r_cut));} end; procedure prepare_phi_numerical; var in_dat : text; name : string; r,ph : array[0..100] of double; i,n,k : integer; x : double; //dummy1 : double; //dummy2 : double; begin write('Enter filename to load : '); // screening function in Angstrom readln(name); assign(in_dat,name); reset(in_dat); r[0] := 0; ph[0] := 1; i := 1; repeat readln(in_dat,r[i],ph[i]); r[i] := r[i]/0.529; inc(i); until eof(in_dat); close(in_dat); n := i-1; rmax := r[n]; writeln('rmax = ',rmax:7:1); for i := 0 to nmax do begin x := i*rmax/nmax; k := -1; repeat inc(k); until (k>n) or ((x>=r[k]) and (x<=r[k+1])); if k>n then halt; phi[i] := ((x-r[k])/(r[k+1]-r[k])*((ph[k+1])-(ph[k])) + (ph[k])); end; end; procedure prepare_phi_scaling(rcut : re; pot_type : byte); var i : integer; x : double; begin rmax := rcut; for i := 0 to nmax do begin x := i*rmax/nmax/a_scr; if pot_type = 1 then begin {moliere screening function} phi[i] := 0.35*exp(-0.3*x)+0.55*exp(-1.2*x)+0.10*exp(-6.0*x); end; if pot_type = 2 then begin {zbl screening function} phi[i] := 0.18175*EXP(-3.1998*x) + 0.50986*EXP(-0.94229*x) + 0.28022*EXP(-0.4029*x) + 0.028171*EXP(-0.20162*x); end; end; end; procedure magic(eps,b:re; var c,r : re); // c = cos(theta/2) eps = E_cm/(Z1 Z2 e^2/a_scr) b = p/a_scr r=r_0/a_scr var rr,v,v1,fr,fr1, ex1,ex2,ex3,ex4, q,roc,sqe,cc,aa,ff, delta :re; begin R:=B; RR:=ln(0.02817/(EPS*B))/0.2016; IF(RR>B) then begin RR:=ln(0.02817/(EPS*RR))/0.2016; IF(RR>B)then R:=RR; end; REPEAT EX1:=0.18175*EXP(-3.1998*R); EX2:=0.50986*EXP(-0.94229*R); EX3:=0.28022*EXP(-0.4029*R); EX4:=0.028171*EXP(-0.20162*R); V:=(EX1+EX2+EX3+EX4)/R; V1:=-(V+3.1998*EX1+0.94229*EX2+0.4029*EX3+0.20162*EX4)/R; FR:=B*B/R+V*R/EPS-R; FR1:=-B*B/(R*R)+(V+V1*R)/EPS-1; Q:=FR/FR1; R:=R-Q; UNTIL (ABS(Q/R)<0.001); ROC:=-2.0*(EPS-V)/V1; SQE:=SQRT(EPS); CC:=(0.01185+SQE)/(0.0068338+SQE); AA:=2.0*EPS*(1.0+(0.80061/SQE))*power(B,CC); FF:=(SQRT(AA*AA+1.0)-AA)*((10.855+EPS)/(16.883+EPS)); DELTA:=(R-B)*AA*FF/(FF+1.0); C:=(B+DELTA+ROC)/(R+ROC); end; var U,W : array[1..16*1024] of double; N_mehler : integer; procedure MEHLER(N : integer); const HALF = 0.5; var W0,Y,Z : double; J : integer; begin N_mehler := N; W0 := PI / (N + N); Y := COS( W0 ); Z := HALF * Y; U[1] := SQRT( HALF + Z ); U[N] := SQRT( HALF - Z ); Z := U[1] * U[N] + U[1] * U[N]; W[1] := W0 * U[N]; W[N] := W0 * U[1]; for J := 2 to (N+1) div 2 do begin U[J] := Y * U[J-1] - Z * U[N-J+2]; U[N-J+1] := Z * U[J-1] + Y * U[N-J+2]; W[J] := W0 * U[N-J+1]; W[N-J+1] := W0 * U[J]; end; end; procedure scattering(var theta : re); var i : integer; q, qa,qb, qx, qq, qz : double; begin get_ro; qz := 0; qx := ro/Ecm; for i := 1 to N_mehler do begin QA := ro/ U[I]; QQ := b* U[I]; QB := qx * U[I] * V_pot(QA)*qa; QQ := (ro + QQ) * (ro - QQ); Q := W[I] / SQRT(QQ - QB); QZ := QZ + Q end; QA := b * QZ; theta := pi-2*qa; end; procedure test_mehler; var theta1,theta2,theta3,c,r0 : re; begin MEHLER(16); b := 1; repeat write('b = ');readln(b); get_theta(theta1); writeln; write('simp = ',theta1*180/pi:8:6); scattering(theta2); write(' mehler = ', theta2*180/pi:8:6); magic(Ecm/(z1*z2/a_scr),b/a_scr,c,r0); theta3 := 2*arccos(c); writeln(' magic(zbl) = ', theta3*180/pi:8:6); until b > 10; end; procedure scan_b(bmax : re; nb : longint; theta_initial,dtheta : double; ntheta : integer ; var cross_ratio : cross_ratio_type); var theta : re; i : longint; db,b_old, b_mean, Rutherford, theta_mean_std, theta_mean_rec, theta_mean, theta_old : re; cross_ratio1:re; itheta : integer; begin MEHLER(16*1024); theta := pi; b := 0; // for itheta := 1 to ntheta do cross_ratio[itheta] := 1e100; <====== for i := 1 to nb do begin theta_old := theta; b_old := b; b := bmax*power(i/nb,3); db := b - b_old; b_mean := 0.5*(b+b_old); // ep := 1e-4; get_theta(theta); // simpson ep := 1e-3; scattering(theta); // mehler // theta := 2*arctan(Z1*Z2/Ecm/b/2); theta_mean_std := 0.5*(theta_old+theta); theta_mean_rec := 1/(0.5*(1/theta_old+1/theta)); theta_mean := 0.6*theta_mean_rec + 0.4*theta_mean_std; Rutherford := sqr(Z1*Z2/4/Ecm) / sqr(sqr(sin(theta_mean/2))); cross_ratio1 := abs(b_mean/sin(theta_mean) * db/(theta_old-theta)) / Rutherford; for itheta := 1 to ntheta do if (theta_initial+(itheta-1)*dtheta <= theta_old) and (theta_initial+(itheta-1)*dtheta >= theta) then begin cross_ratio[itheta] := cross_ratio1; // <==== end; end; end; procedure get_cross_ratio(Z1_in,m1_in,Z2_in,m2_in,q_in,E_in,b_max_in : double; nb_in: integer; ep_in : double; theta_initial,dtheta : double; ntheta : integer; var cross_ratio : cross_ratio_type); var b_max, r_max, m1,m2,E,q : double; begin z1 := Z1_in {6} {79} {1, 2, 14, 92}; // H, He, Si, U q := q_in {3} {46}; z2 := z2_in {13}; // Si 14, Hf 72 m1 := m1_in {12} {197} {1, 4, 28, 238}; m2 := m2_in {27}; // Si 28, Hf 178 ep := ep_in {1e-5}; // prepare_phi_numerical; b_max := b_max_in {1}; r_max := 10*b_max; a_scr := 0.8854/sqrt(power(Z1-q,2/3)+power(Z2,2/3)); prepare_phi_scaling(r_max,1); E := E_in*m1; // 25, 100, 1000 5000 keV/u E := E*1000/27.2; Ecm := E*m2/(m1+m2); Tmax := 4*m1*m2/sqr(m1+m2)*E; vmax := sqrt(2*Tmax/(m2*1836)); // test_mehler; scan_b(b_max,nb_in,theta_initial, dtheta,ntheta,cross_ratio); end; end.