c......Example 5(b): Modified Cartesian to spherical polar coordinates c..... convertor: Part 5 c23456 program ex5b c implicit none c double precision cart(3),polar(3) integer ichoice c c..... Declare the two functions c double precision degtorad,radtodeg c c..... Read in whether to convert Cartesian to spherical polar c..... or to convert spherical polar to Cartesian c c....NOTE: throughout this program I have removed the print statements c.... needed if the program is run interactively using the screen c.... and keyboard c read(5,*) ichoice c c..... Beginning of if loop c c..... Perform Cartesian to spherical polar conversion c if(ichoice.eq.1) then c c..... Read in the three Cartesian coordinates from the file c read(5,*) cart(1),cart(2),cart(3) c c..... conversion c call ctosp(cart,polar) c c..... Remember that the trigonometric units are by default radians c..... To convert to degrees you need a value for pi, c..... Recall: 2*pi radians = 360 degrees c polar(2) = radtodeg(polar(2)) polar(3) = radtodeg(polar(3)) c c..... Output the result to the keyboard c..... The angular results are output in degrees c write(6,*) 'Spherical polar coordinates are (angles in degrees)' write(6,*) ' r = ',polar(1) write(6,*) ' theta = ',polar(2) write(6,*) ' phi = ',polar(3) c c..... Spherical polar to cartesian conversion c else if(ichoice.eq.2) then c c..... Read in the three spherical polar coordiantes from a file c read(5,*) polar(1),polar(2),polar(3) c c..... Convert angles in degrees to angles in radians c polar(2) = degtorad(polar(2)) polar(3) = degtorad(polar(3)) c c..... Call the subroutine to do the spherical polar to Cartesian c..... conversion c call sptoc(polar,cart) c c..... Output the results c write(6,*) 'Cartesian coordiantes are:' write(6,*) 'x = ',cart(1) write(6,*) 'y = ',cart(2) write(6,*) 'z = ',cart(3) c c..... If 1 or 2 is not entered print a warning message else write(6,*) 'You MUST enter 1 or 2!' c c..... End of the if-else block c end if c end c c..... Subroutine to calculate from Cartesian to spherical polar c..... coordinates c23456 subroutine ctosp(x,y) c double precision x(3),y(3) c c.....NOTE: There are a variety of ways of converting from Cartesian to c..... spherical polar coordinates. The method illustratd below is c..... just one. c c..... Calculate the r vector c..... Assign the r vector to polar(1) c y(1) = sqrt(x(1)**2+x(2)**2+x(3)**2) c c..... Calculate the value of theta using the fact that z = r*cos(theta) c..... Assign theta to polar(2) c y(2) = acos(x(3)/y(1)) c c..... Calculate the value of phi using the fact that c..... x = r*sin(theta)*cos(phi) c y(3) = acos(x(1)/(y(1)*sin(y(2)))) c c..... Return to main program c return c end c c..... Subroutine to convert from spherical polar coordinates to Cartesian c..... coordinates c subroutine sptoc(y,x) c double precision y(3),x(3) c c..... Perform the simple conversion using the well known formulae c..... x = r*sin(theta)*cos(phi) c..... y = r*sin(theta)*sin(phi) c..... z = r*cos(theta) c x(1) = y(1)*sin(y(2))*cos(y(3)) x(2) = y(1)*sin(y(2))*sin(y(3)) x(3) = y(1)*cos(y(2)) c c..... Retun to the main program return end c c.... Function to convert degrees to radians double precision function degtorad(x) c implicit none double precision x,pi c c.... Define pi as a parameter parameter(pi = 3.14159265358979d0) c degtorad = x *pi/180.0d0 c return end c c.... Function to convert radians to degrees double precision function radtodeg(x) c implicit none double precision x,pi c c.... Define pi as a parameter parameter(pi = 3.14159265358979d0) c radtodeg = x *180.0d0/pi c return end
This page maintained by alex.brown@ualberta.ca of the Department of Chemistry, University of Alberta
Last updated August 8, 2003.