c......Example 5(a): Modified Cartesian to spherical polar coordinates
c.....               convertor: Part 4
c23456
      program ex5
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..... Open the input and output files
c
      open(unit = 7,file = 'input.dat')
      open(unit = 8,file = 'output.dat')
c
c..... some compilers require the STATUS of the file to be declared
c..... and you will have to use
c      open(unit = 7,file = 'input.dat',status='old')
c      open(unit = 8,file = 'output.dat',status='new')
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(7,*) 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(7,*) 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(8,*) 'Spherical polar coordinates are (angles in degrees)'
       write(8,*) ' r = ',polar(1)
       write(8,*) ' theta = ',polar(2)
       write(8,*) ' 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(7,*) 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(8,*) 'Cartesian coordiantes are:'
         write(8,*) 'x = ',cart(1)
         write(8,*) 'y = ',cart(2)
         write(8,*) 'z = ',cart(3)
c
c..... If 1 or 2 is not entered print a warning message
       else
         write(8,*) '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