ENTER A or 0 for Schwarzschild radius for our sun
0.000000000000000E+000
MASS OF BODY - KILOGRAMS
1.990000000000000E+030
RADIUS OF THE SUN - METERS
695950000.000000
TOTAL RAY DEFLECTION IS 1.75074437282132 arc seconds official
WIDTH OF LIGHT RAY - METERS
655000.000000000
TIME INCREMENT - SECONDS
1.000000000000000E-002
1.28462126313333 ARC SEC, DIST = 1 R
1.59054184508388 ARC SEC, DIST = 2 R
1.67527376422873 ARC SEC, DIST = 3 R
1.70827290038373 ARC SEC, DIST = 4 R
1.72421867880825 ARC SEC, DIST = 5 R
1.73306790110479 ARC SEC, DIST = 6 R
1.73846825118068 ARC SEC, DIST = 7 R
1.74200061229695 ARC SEC, DIST = 8 R
1.74443500433983 ARC SEC, DIST = 9 R
1.74618346603163 ARC SEC, DIST = 10 R
1.74747904551263 ARC SEC, DIST = 11 R
1.74845208156234 ARC SEC, DIST = 12 R
1.74921369249508 ARC SEC, DIST = 13 R
1.74982033842731 ARC SEC, DIST = 14 R
1.75031166156846 ARC SEC, DIST = 15 R
1.75072850540704 ARC SEC, DIST = 16 R
TOTAL RAY DEFLECTION IS 1.75072850540704 ARC SEC
* TOTAL RAY DEFLECTION IS 1.75074437282132 arc seconds official for standard model
* in 12/5/2023 I get defl 1.75072850540704 arc seconds using 2GM/c^2 M=sun mass
* bending of light around the sun using epi = epi0*(1+A/R)**2 bending Dec 5, 2023
* Ray starts at the sun and leaves the sun.
* The angle of departure is doubled to get the total angle for both arrival and departure angles.
* X0 Y0 IS ON THE SURFACE OF THE SUN, X1 Y1 IS FARTHER OFF THE SURFACE
IMPLICIT DOUBLE PRECISION (A-H,L-Z)
* the 4000 here results in a distance away from the Sun of 8 diameters
PRINT *,'ENTER A or 0 for Schwarzschild radius for our sun'
READ *,A !
PRINT *,A
G=6.67D-11
PI=3.14159265D0
Co=2.998D8
EPI0=8.85418783D-12
MU=4.*3.14159265D-7
PRINT *,'MASS OF BODY - KILOGRAMS'
M=1.99D30
* READ *,M
IF(A.LE.0.) A=2*G*M/Co**2 ! Schwarzschild radius formula
WRITE(*,*) M
PRINT *,'RADIUS OF THE SUN - METERS'
Radius=695950.D3
* READ *,Radius
WRITE(*,*) Radius
* below is a check for the bending around the sun using an official calculation
* TOTAL RAY DEFLECTION IS 1.75074437282132 arc seconds official due to general relativity
* REL13 RAY DEFLECTION IS 1.75072850540704 arc seconds using 2GM/c^2 M=sun mass
ANG=(4.*G*M/Radius/Co**2)*360.*3600./2./PI ! official answer to the angle
PRINT *,'TOTAL RAY DEFLECTION IS',ANG,' arc seconds official'
* this width of the beam produces the same bending as the exact formula
W=655000.D0
3 PRINT *,'WIDTH OF LIGHT RAY - METERS'
* READ *,W
PRINT *,W
* ~750,000 km radius in ~2.5 sec. if there are 250 steps each radius distance time, .01 sec each step then 4000 of these
* is a total time of 40 sec or 4000/250 = 16 radius's or 8 sun diameters the wave travels
TI = .01D0
PRINT *,'TIME INCREMENT - SECONDS'
* READ *,TI
PRINT *,TI
* the angle must be small and the ray is at Y0 bottom of the ray which is traveling along X to the right
Y0=Radius
Y1=Radius+W
* .01 sec per interval produces a TT total time of 40 sec or 16 times the radius of the Sun
TT=0.D0
X0=0.D0
X1=0.D0
DO 1 I=1,4000
* ray travels in X direction only and assume change in Y is negligible (this is an approximation, good only for very small angles)
R0=DSQRT(X0*X0 + Y0*Y0)
R1=DSQRT(X1*X1 + Y1*Y1)
* formula EPI = EPI0*(1+A/R)**2 is tested
* deviation in distance with respect to 0
EPI = EPI0*(1.D0+A/R0)**2
X00 = TI/(DSQRT(MU*EPI))
EPI = EPI0*(1.D0+A/R1)**2
X11 = TI/(DSQRT(MU*EPI))
* remember the distance traveled by the two rays as we move away from the sun
X0=X0+X00
X1=X1+X11
* XD is the difference in the small deviations, note that subtraction above is reflected in the reverse of the subtraction below
XD=X1-X0
IF(I/250*250.EQ.I)
& PRINT *,2.D0*360.D0*3600.D0*XD/(2.*PI*W),
& ' ARC SEC, DIST =',I/250,' R'
1 CONTINUE
* the factor of 2 is to include both approach and leaving of the wave, i.e. double the leaving angle
ANG=2.D0*360.D0*3600.D0*XD/(2.*PI*W)
PRINT *,'TOTAL RAY DEFLECTION IS',ANG,' ARC SEC'
* IF(W.GT.10.) GOTO 3
STOP
END
An executable for this program is posted here:
https://egpreston.com/REL13.exe and it's not a virus.
If successfully downloaded put it in a folder where you can open a command prompt window
and type dir to make sure you can see the REL13.exe file and then type REL13.exe and it
will run the program and you can see the results are the same as at the top of this file
if you enter 0 for the value of A. Or you can copy the program and use your own compiler.