! Physics 222: Fortran/C++ Programing for Scientist
! Written by Adam Schendel started Feb 15th, 2007

!  This program will calculate where a light beam will strike a object and how the light beam will 
!  react in response to the interaction.  New objects can also be created in this program.  

!jc It would be nice to have a little longer introduction to the program here.


! Fortran 90/95 version
!Complile with gfortran -o light light.f90 -lpgplot

MODULE light_mods
!jc  A little explanation of the variables would be nice - purpose, units, etc.
   TYPE::prism3                 !This module is used to hold info about three sided objects
      REAL::refrac_index, side(3,4), vertex(3,2), scale, x_offset, y_offset, theta_offset
      LOGICAL::invert
    END TYPE

    TYPE::prism4                 !This module is used to hold info about four sided objects 
      REAL::refrac_index, side(4,4), vertex(4,2), scale, x_offset, y_offset, theta_offset
      LOGICAL::invert
    END TYPE

    TYPE::cross                  !This module is used to hold info about intersection points
      REAL::refrac_index, side(4), intersect(2)
    END TYPE
END MODULE       

!jc  Nice use of modules and structures.


PROGRAM light
USE light_mods
IMPLICIT NONE

CHARACTER(len=22)::prism_type
CHARACTER(len=20):: option, title, charac_i
INTEGER::prism3_num, prism4_num, new_prism, invert_angle, i, j, k, insec_num, ierror, refrac_num
INTEGER(KIND=2):: prism3_angle(3), prism4_angle(4), ier, pgbeg
REAL:: length3(3), length4(4), vertex(5,2), all_intersect(2), beam(1000,4), true_intersect(1001,2)
REAL:: beam_angle(1000), refractive(2), slope, table(4)
LOGICAL::touch, bound
TYPE (prism3):: prism3_(20)
TYPE (prism4):: prism4_(20)
TYPE (cross)::test_intersect(200)
TYPE (prism4)::boundary


    10 FORMAT (2F6.2) !Default write vertex format used to write vertex points
    WRITE(*,*) '	This program will calculate where a light beam will strike'
    WRITE(*,*) 'a object and how the light beam will react in response to the'
    WRITE(*,*) 'interaction.  New objects can also be created in this program.'
    WRITE(*,*) ''
    WRITE(*,*) 'Would you like to create a new prism? [Y/N]'
    READ(*,*) option  
    !jc  This odd indentation - indent inside the if block, not the if itself
      IF ((option == 'Y') .OR. (option == 'y') .OR. (option == 'Yes') .OR. (option == 'yes')) THEN
        DO i=1,3        !This if statement lets the user create 3-sided prisms and gives three chances to get it right

          !This promps for the angles of the prism
          WRITE(*,*) ''
          WRITE(*,*) 'ALL ANGLES IN THIS PRISM MUST ADD UP TO 180 DEGREES!!'
          WRITE(*,*) 'What is the angle between side 3 and side 1?'
          READ(*,*) prism3_angle(1)
          WRITE(*,*) 'What is the angle between side 1 and side 2?'
          READ(*,*) prism3_angle(2)
          WRITE(*,*) 'What is the angle between side 2 and side 3?'
          READ(*,*) prism3_angle(3)

          !This promps for the position of each vertex which are then used to calculate the sides.
          WRITE(*,*) 'What is the position of each of the three vertex in cm. Example X,Y'
          WRITE(*,*) 'Vertex 1'
          READ(*,*) prism3_(1)%vertex(1,1:2)
          WRITE(*,*) 'Vertex 2'
          READ(*,*) prism3_(1)%vertex(2,1:2)
          WRITE(*,*) 'Vertex 3'
          READ(*,*) prism3_(1)%vertex(3,1:2)
          CALL sides3(prism3_(1)%vertex(1:3,1:2), prism3_(1)%side(1:3,1:4))      !This creates the equations for the sides of the prism

!jc  Error checking on input would be nice, to make this more robust.

            WRITE(*,*) 'Was this prism created correctly? [Y/N]'
            READ(*,*) option  
            IF ((option == 'Y') .OR. (option == 'y') .OR. (option == 'Yes') .OR. (option == 'yes')) THEN
              !Creating a three side prism file
              WRITE(*,*) 'What is the name of the file that you wish to save this prism as?'
              WRITE(*,*) 'An exmaple of a good file name is prism3_0306090.f90 where 0306090 refer to'
              WRITE(*,*) 'the angles of the prism in asending order.'
              READ(*,*) prism_type
              OPEN(UNIT=11, FIlE='prisms/'//prism_type)
              WRITE(11,*) prism3_(1)%vertex(1,1:2)
              WRITE(11,*) prism3_(1)%vertex(2,1:2)
              WRITE(11,*) prism3_(1)%vertex(3,1:2)
              CLOSE(11)
            
              !Creating the length of the sides of the prism
              length3(1)=SQRT(((prism3_(1)%vertex(1,1))-(prism3_(1)%vertex(2,1)))**2+&
              ((prism3_(1)%vertex(1,2))-(prism3_(1)%vertex(2,2)))**2)
              length3(2)=SQRT(((prism3_(1)%vertex(2,1))-(prism3_(1)%vertex(3,1)))**2+&
              ((prism3_(1)%vertex(2,2))-(prism3_(1)%vertex(3,2)))**2)
              length3(3)=SQRT(((prism3_(1)%vertex(3,1))-(prism3_(1)%vertex(1,1)))**2+&
              ((prism3_(1)%vertex(3,2))-(prism3_(1)%vertex(1,2)))**2)
                        
              !This creates a entry in the file prism3_option.f90
              OPEN(UNIT=12, FILE='prisms/prism3_option.f90', STATUS='OLD', ACCESS='APPEND')
              WRITE(12,*) prism_type 
              WRITE(12,100) prism3_angle(1:3), length3(1:3)
              100 FORMAT ('                      ', 3I6, '    ', 3F6.2)
              CLOSE(12)
              EXIT
            END IF
        END DO
      END IF

    WRITE(*,*) ''
    WRITE(*,*) 'Would you like to perform a run? [Y/N]'
    option='hi'
    READ(*,*) option  
    IF ((option == 'Y') .OR. (option == 'y') .OR. (option == 'Yes') .OR. (option == 'yes')) THEN
 
 !jc  When you do the same sort of thing over and over, it is nice to abstract
 !jc  it into a subroutine.  You could have a subroutine or function for y/n
 !jc  questions.  Send in a string for the questions to ask, and get back
 !jc  the answer.

      !This section defines the size of the table used that holds the prisms
      WRITE(*,*) ''
      WRITE(*,*) 'Relative to the point 0,0 where are the four sides of the table located?'
      WRITE(*,*) 'Table dimensions should be rounded to the nearest centimeter.'
      WRITE(*,*) 'Where is the left side of the table?'
      READ(*,*) table(1)
      WRITE(*,*) 'Where is the right side of the table?'
      READ(*,*) table(2)
      WRITE(*,*) 'Where is the far side of the table?'
      READ(*,*) table(3)
      WRITE(*,*) 'Where is the near side of the table?'
      READ(*,*) table(4)

      boundary%vertex(1,1)=table(1)
      boundary%vertex(1,2)=table(4)
      boundary%vertex(2,1)=table(2)
      boundary%vertex(2,2)=table(4)
      boundary%vertex(3,1)=table(2)
      boundary%vertex(3,2)=table(3)
      boundary%vertex(4,1)=table(1)
      boundary%vertex(4,2)=table(3)

      CALL sides4(boundary%vertex(1:4,1:2), boundary%side(1:4,1:4))   
        
      WRITE(*,*) 'What is the refractive index of the empty space on the table?'
      READ(*,*) refractive(1)
      
      
      ! This portion of the program performs the individual runs    
      WRITE(*,*) ''
      WRITE(*,*) 'How many three sided prism do you want to use in this run?'
      READ(*,*)prism3_num   !This section sets the number of prisms and then set the type and position
      DO i=1,prism3_num
        CALL SYSTEM('cat ' // 'prisms/prism3_option.f90')
                  !jc  You could just combine the two strings and get rid of //
        WRITE(*,*) ''
        WRITE(*,*) 'What type of three sided prism would you like to use as'
        WRITE(*,102) i
        102 FORMAT ('prism3_ ', I1,'.  Example "prism3_0306090.f90".')
        READ(*,*) prism_type
        OPEN(UNIT=10, FILE='prisms/'//prism_type)
        DO j=1,3
          READ(10,*) prism3_(i)%vertex(j,1:2)
        END DO
        CLOSE(UNIT=10)
        WRITE(*,*) 'What is the scale of the prism?  Example scale 4 on a 30x60x90 will' 
        WRITE(*,*) 'make a prism with sides of length 4cm, 2.32cm and 4.6cm.'
        READ(*,*) prism3_(i)%scale
        DO j=2,3
          DO k=1,2
            prism3_(i)%vertex(j,k)=prism3_(i)%vertex(1,k)+(prism3_(i)%scale)*((prism3_(i)%vertex(j,k))-(prism3_(i)%vertex(1,k)))
          END DO
        END DO
        
        !This section inverts the prism about the y-axis
        WRITE(*,*) 'Would you like to invert the prism about the y-axis? [Y/N]'
        prism3_(i)%invert=.false.
        READ(*,*) option
        IF ((option == 'Y') .OR. (option == 'y') .OR. (option == 'Yes') .OR. (option == 'yes')) THEN
          prism3_(i)%invert=.true.
          DO j=1,3
            prism3_(i)%vertex(j,1)=(-1)*(prism3_(i)%vertex(j,1))

          END DO                       !This section find the intersection and then finds the new beam 
          vertex(1,1)=prism3_(i)%vertex(2,1)+ABS((prism3_(i)%vertex(3,1))-(prism3_(i)%vertex(1,1)))
          vertex(2,2)=prism3_(i)%vertex(2,2)
          vertex(2,1)=prism3_(i)%vertex(1,1)+ABS((prism3_(i)%vertex(3,1))-(prism3_(i)%vertex(1,1)))
          vertex(2,2)=prism3_(i)%vertex(1,2)
          vertex(3,1)=prism3_(i)%vertex(3,1)+ABS((prism3_(i)%vertex(3,1))-(prism3_(i)%vertex(1,1)))
          vertex(3,2)=prism3_(i)%vertex(3,2)
          !jc  You could have shortened this by using array syntax:
          !jc  vertex(:,2)=prism_(i)%vertex(:2)
          !jc  etc.
          DO j=1,3
            WRITE(*,10) vertex(j,1:2)
            prism3_(i)%vertex(j,1:2)=vertex(j,1:2)
          END DO
          
        END IF        
           
        !This section rotates the prism about the first vertex
        WRITE(*,*) 'What is the theta offset of the prism from the x axis in degrees?'
        READ(*,*) prism3_(i)%theta_offset
        vertex(1:3,1:2)=prism3_(i)%vertex(1:3,1:2)
        
       !jc  becareful using integers in floating point calculations 
        prism3_(i)%vertex(1,1:2)=prism3_(i)%vertex(1,1:2)
        prism3_(i)%vertex(2,1)=SQRT((vertex(2,1)-vertex(1,1))**2+(vertex(2,2)-vertex(1,2))**2)&
          *COS((prism3_(i)%theta_offset*2*(3.14159265)/360+&
          ATAN(((vertex(2,2))-(vertex(1,2)))/((vertex(2,1))-(vertex(1,1))))))
        prism3_(i)%vertex(2,2)=SQRT((vertex(2,1)-vertex(1,1))**2+(vertex(2,2)-vertex(1,2))**2)&
          *SIN((prism3_(i)%theta_offset*2*(3.14159265)/360+&
          ATAN(((vertex(2,2))-(vertex(1,2)))/((vertex(2,1))-(vertex(1,1))))))
        prism3_(i)%vertex(3,1)=SQRT((vertex(3,1)-vertex(1,1))**2+(vertex(3,2)-vertex(1,2))**2)&
          *COS((prism3_(i)%theta_offset*2*(3.14159265)/360+&
          ATAN(((vertex(3,2))-(vertex(1,2)))/((vertex(3,1))-(vertex(1,1))))))
        prism3_(i)%vertex(3,2)=SQRT((vertex(3,1)-vertex(1,1))**2+(vertex(3,2)-vertex(1,2))**2)&
          *SIN((prism3_(i)%theta_offset*2*(3.14159265)/360+&
          ATAN(((vertex(3,2))-(vertex(1,2)))/((vertex(3,1))-(vertex(1,1))))))
          
        WRITE(*,10) prism3_(i)%vertex(1,1:2)
        WRITE(*,10) prism3_(i)%vertex(2,1:2)
        WRITE(*,10) prism3_(i)%vertex(3,1:2)

        !This section locates the position of the firts vertex
        WRITE(*,*) 'What is the position of the first vertex.'
        WRITE(*,*) 'Example "X,Y".'
        READ(*,*) prism3_(i)%x_offset, prism3_(i)%y_offset
        DO j=1,3
          prism3_(i)%vertex(j,1)=(prism3_(i)%vertex(j,1))+(prism3_(i)%x_offset)
          prism3_(i)%vertex(j,2)=(prism3_(i)%vertex(j,2))+(prism3_(i)%y_offset)
        END DO
        WRITE(*,10) prism3_(i)%vertex(1,1:2)
        WRITE(*,10) prism3_(i)%vertex(2,1:2)
        WRITE(*,10) prism3_(i)%vertex(3,1:2)
        
        !This section set the refractive index of the prism        
        WRITE(*,103) i
        103 FORMAT (' What is the refractive index of prism3_', I2)
        READ(*,*) prism3_(i)%refrac_index
        CALL sides3(prism3_(i)%vertex(1:3,1:2), prism3_(i)%side(1:3,1:4)) 
      
      END DO

      DO
        WRITE(*,*) 'Would you like perform a run with these prisms? [Y/N]'
        READ(*,*) option  
        IF ((option == 'Y') .OR. (option == 'y') .OR. (option == 'Yes') .OR. (option == 'yes')) THEN
    
        !This section defines the equation of the first beam section.
        WRITE(*,*) 'What is the origin of the light beam?  Example X,Y'
        READ(*,*) true_intersect(1,1:2)
        WRITE(*,*) 'What is the angle in degrees that the light beam makes'
        WRITE(*,*) 'with the positive x-axis? Example 45'
        READ(*,*) beam_angle(1)

        ! Initialize plot
        ier = pgbeg(0, '/xserve',1,1)
        ! set scale for plot
        CALL pgenv(table(1), table(2), table(4), table(3), 1, 1)
        !set labels
        CALL pglab ('', '', 'Light Table Simulation')
        !origin of beam
        CALL pgpt1(true_intersect(1,1),true_intersect(1,2),0565)
        insec_num=0
        DO 
          OPEN(UNIT=14, FILE='prisms/test_intersects.f90', IOSTAT=ierror)
          insec_num=insec_num+1
        
          !This section defins the equation of the beam that will be tested
          DO j=1,10
            IF (beam_angle(insec_num)==(90*(2*j-1))) THEN
              beam(insec_num,1)=0
              beam(insec_num,2)=1
              beam(insec_num,3)=1
              beam(insec_num,4)=true_intersect(insec_num,1)
              EXIT
            ELSE IF (beam_angle(insec_num)==(90*(2*j-2))) THEN
              beam(insec_num,1)=1
              beam(insec_num,2)=0
              beam(insec_num,3)=0
              beam(insec_num,4)=true_intersect(insec_num,2)
              EXIT
            ELSE
              beam(insec_num,1)=1
              beam(insec_num,2)=(TAN(beam_angle(insec_num)*2*(3.14159265)/360))
              beam(insec_num,3)=1
              beam(insec_num,4)=((true_intersect(insec_num,2))-(beam(insec_num,2))*(true_intersect(insec_num,1)))
            END IF
          END DO

          DO i=1, prism3_num    
            vertex(1,1:2)=prism3_(i)%vertex(1,1:2)
            vertex(2,1:2)=prism3_(i)%vertex(2,1:2)
            vertex(3,1:2)=prism3_(i)%vertex(3,1:2)
            vertex(4,1:2)=prism3_(i)%vertex(1,1:2)

            !This plots the boundary of the prism3
            CALL PGLINE(4,vertex(1:4,1),vertex(1:4,2))
            !Position of prism lables
            WRITE(charac_i, '(i1)') i
            title='Prism3_' // charac_i
            CALL PGPTXT(prism3_(i)%vertex(1,1),prism3_(i)%vertex(1,2),0.0,0.0,title)

           !This section finds all possible intersections
            DO j=1,3
              CALL intersect(beam(insec_num,1:4), prism3_(i)%side(j,1:4), vertex(j:j+1,1:2), all_intersect(1:2), touch)
              WRITE(*,*) touch
              touch=.FALSE.

              IF  (((vertex(j,1)-0.1<all_intersect(1)).AND.(vertex(j+1,1)+0.1>all_intersect(1))).OR. &
              ((vertex(j,1)+0.1>all_intersect(1)).AND.(vertex(j+1,1)-0.1<all_intersect(1)))) THEN
                IF (((vertex(j,2)+0.1>all_intersect(2)).AND.(vertex(j+1,2)-0.1<all_intersect(2))).OR. &
                  ((vertex(j,2)-0.1<all_intersect(2)).AND.(vertex(j+1,2)+0.1>all_intersect(2)))) THEN

                  touch=.TRUE.

                  WRITE(14,*) prism3_(i)%refrac_index, prism3_(i)%side(j,1:4), all_intersect(1:2)
                END IF 
              END IF
              WRITE(*,*) ''

            END DO
          END DO

          !This section finds the intersection along the boundary
          vertex(1:4,1:2)=boundary%vertex(1:4,1:2)
          vertex(5,1:2)=boundary%vertex(1,1:2)        
          DO j=1,4
            CALL intersect(beam(insec_num,1:4), boundary%side(j,1:4), vertex(j:j+1,1:2), all_intersect(1:2), bound)
            WRITE(14,*) boundary%refrac_index, boundary%side(j,1:4), all_intersect(1:2)
            WRITE(*,*) boundary%side(j,1:4)
          END DO
        
          CLOSE(14)
      
          !This section finds the true intersection
          CALL true_intersection(true_intersect(insec_num,1:2), beam_angle(insec_num), true_intersect(insec_num+1,1:2),&
            refractive(2), slope)

          !This section determins how the beam will act when it hits the different surfaces
          IF (MOD(refrac_num,2)/=0) THEN
            CALL new_beam(beam_angle(insec_num), slope, refractive(2), refractive(1), refrac_num,&
              beam_angle(insec_num+1), refrac_num)
          ELSE IF (MOD(refrac_num,2)==0) THEN
            CALL new_beam(beam_angle(insec_num), slope, refractive(1), refractive(2), refrac_num,&
              beam_angle(insec_num+1), refrac_num)
          END IF

          IF  ((true_intersect(insec_num,1)<=table(1)) .OR. (true_intersect(insec_num,1)>=table(2)) .OR. &
          (true_intersect(insec_num,2)>=table(3)) .OR. (true_intersect(insec_num,2)<=table(4))) THEN
            WRITE(*,*) 'The beam has left the table.'
            EXIT
          END IF          
          WRITE(*,*) '=========================================='
        END DO

        DO i=1,insec_num
          WRITE(*,*) beam_angle(i)
          WRITE(*,*) true_intersect(i,1:2)
          WRITE(*,*) beam(i,1:4)
        END DO
        WRITE(*,*) insec_num

        !Plot the beam path
        CALL PGLINE(insec_num,true_intersect(1:insec_num,1),true_intersect(1:insec_num,2))

        !close the plot
        CALL pgend 
        ELSE
          EXIT
        END IF

        !This writes a list prisms refractive indexes
        WRITE(*,*)'=============================================='
        WRITE(*,*)' Prism          |     Refractive Index'
        WRITE(*,*)'=============================================='
        DO j=1, prism3_num
          WRITE(*,300) j, prism3_(j)%refrac_index
          300 FORMAT ('  Prism3_', I2, '     |     ', F5.2)
        END DO
        WRITE(*,*)'=============================================='
        
      END DO
    END IF

STOP
END PROGRAM


SUBROUTINE new_beam(beam_angle1, slope, refractive1, refractive2, refrac_num1, beam_angle2, refrac_num2)
IMPLICIT NONE

REAL, INTENT(IN)::beam_angle1, slope, refractive1, refractive2
REAL, INTENT(OUT)::beam_angle2
INTEGER, INTENT(IN)::refrac_num1
INTEGER, INTENT(OUT)::refrac_num2
REAL::size_angle, test_angle1, test_angle2, reflect
INTEGER::i
      
    WRITE(*,*) beam_angle1
    WRITE(*,*) slope
    WRITE(*,*)'-----------'

    WRITE(*,*) refractive1
    WRITE(*,*) refractive2

    size_angle=beam_angle1          !This reduces the angle of the beam to an equivulent value between 0 and 360
    IF (size_angle>=360) THEN
      size_angle=size_angle-360
    END IF

    !This reduces the angle between the beam and the surface of intersection to an equivulent value between 0 and 90
    test_angle1=ABS(90-ABS(((ATAN(slope))*360/(2*(3.14159265)))-size_angle))
    IF (test_angle1 >= 360) THEN
      test_angle1=(test_angle1-360) 
    ELSE IF (test_angle1 >= 270) THEN
      test_angle1=(test_angle1-270)
    ELSE IF (test_angle1 >= 180) THEN
      test_angle1=(test_angle1-180)
    ELSE IF (test_angle1 >= 90) THEN
      test_angle1=(test_angle1-90)
    END IF

    WRITE(*,*) test_angle1

    reflect=(refractive1/refractive2)*(SIN(test_angle1*2*(3.14159265)/360))   !This specifies if the beam will reflect or refract
    IF ((reflect>1.0) .OR. (reflect<-1.0) .OR. (refractive2==(-1.0))) THEN

      !This section governs how the beam reactes is the beam is to reflect
      IF ((((size_angle>0) .AND. (size_angle<=90))) .AND. (slope>10000)) THEN
        beam_angle2=size_angle+(90-test_angle1)*2
      ELSE IF ((((size_angle>0) .AND. (size_angle<=90))) .AND. (slope>0.0)) THEN
        beam_angle2=size_angle-(90-test_angle1)*2
       ELSE IF ((((size_angle>0) .AND. (size_angle<=90))) .AND. (slope<=0.0)) THEN
        beam_angle2=size_angle+(90-test_angle1)*2
      ELSE IF ((((size_angle>90) .AND. (size_angle<=180))) .AND. (slope>0.0)) THEN
        beam_angle2=size_angle-(90-test_angle1)*2
      ELSE IF ((((size_angle>90) .AND. (size_angle<=180))) .AND. (slope<=0.0)) THEN
        beam_angle2=size_angle+(90-test_angle1)*2
      ELSE IF ((((size_angle>180) .AND. (size_angle<=270))) .AND. (slope>0.0)) THEN
        beam_angle2=size_angle+(90-test_angle1)*2
      ELSE IF ((((size_angle>180) .AND. (size_angle<=270))) .AND. (slope<=0.0)) THEN
        beam_angle2=size_angle-(90-test_angle1)*2
      ELSE IF ((((size_angle>270) .AND. (size_angle<=360))) .AND. (slope>0.0)) THEN
        beam_angle2=size_angle+(90-test_angle1)*2
      ELSE IF ((((size_angle>270) .AND. (size_angle<=360))) .AND. (slope<=0.0)) THEN
        beam_angle2=size_angle+(90-test_angle1)*2
      END IF
WRITE(*,*) 'reflect'

!jc  It seems like there must be a way to simplify this logic.
    ELSE
      !This section governs how the beam reactes is the beam is to reflect
      IF ((((size_angle>0) .AND. (size_angle<=90))) .AND. (slope<-1000)) THEN
        beam_angle2=size_angle-test_angle1+((ASIN((refractive1/refractive2)*&
        (SIN(test_angle1*2*(3.14159265)/360))))*360/(2*3.14159265))
      ELSE IF ((((size_angle>0) .AND. (size_angle<=90))) .AND. (slope>0)) THEN
        beam_angle2=size_angle-test_angle1+((ASIN((refractive1/refractive2)*&
        (SIN(test_angle1*2*(3.14159265)/360))))*360/(2*3.14159265))
      ELSE IF ((((size_angle>0) .AND. (size_angle<=90))) .AND. (slope<=0)) THEN
        beam_angle2=size_angle+test_angle1-((ASIN((refractive1/refractive2)*&
        (SIN(test_angle1*2*(3.14159265)/360))))*360/(2*3.14159265))
      ELSE IF ((((size_angle>90) .AND. (size_angle<=180))) .AND. (slope<-1000)) THEN
        beam_angle2=size_angle-test_angle1+((ASIN((refractive1/refractive2)*&
        (SIN(test_angle1*2*(3.14159265)/360))))*360/(2*3.14159265))
      ELSE IF ((((size_angle>90) .AND. (size_angle<=180))) .AND. (slope>0)) THEN
        beam_angle2 =size_angle-test_angle1+((ASIN((refractive1/refractive2)*&
        (SIN(test_angle1*2*(3.14159265)/360))))*360/(2*3.14159265))
      ELSE IF ((((size_angle>90) .AND. (size_angle<=180))) .AND. (slope<=0)) THEN
        beam_angle2 =size_angle+test_angle1-((ASIN((refractive1/refractive2)*&
        (SIN(test_angle1*2*(3.14159265)/360))))*360/(2*3.14159265))
      ELSE IF ((((size_angle>180) .AND. (size_angle<=270))) .AND. (slope<-1000)) THEN
        beam_angle2=size_angle-test_angle1+((ASIN((refractive1/refractive2)*&
        (SIN(test_angle1*2*(3.14159265)/360))))*360/(2*3.14159265))
      ELSE IF ((((size_angle>180) .AND. (size_angle<=270))) .AND. (slope>0)) THEN
        beam_angle2 =size_angle+test_angle1-((ASIN((refractive1/refractive2)*&
        (SIN(test_angle1*2*(3.14159265)/360))))*360/(2*3.14159265))
      ELSE IF ((((size_angle>180) .AND. (size_angle<=270))) .AND. (slope<=0)) THEN
        beam_angle2 =size_angle+test_angle1+((ASIN((refractive1/refractive2)*&
        (SIN(test_angle1*2*(3.14159265)/360))))*360/(2*3.14159265))
      ELSE IF ((((size_angle>270) .AND. (size_angle<=360))) .AND. (slope<-1000)) THEN
        beam_angle2=size_angle-test_angle1+((ASIN((refractive1/refractive2)*&
        (SIN(test_angle1*2*(3.14159265)/360))))*360/(2*3.14159265))
      ELSE IF ((((size_angle>270) .AND. (size_angle<=360))) .AND. (slope>0)) THEN
        beam_angle2 =size_angle-test_angle1+((ASIN((refractive1/refractive2)*&
        (SIN(test_angle1*2*(3.14159265)/360))))*360/(2*3.14159265))
      ELSE IF ((((size_angle>270) .AND. (size_angle<=360))) .AND. (slope<=0)) THEN
        beam_angle2 =size_angle-test_angle1+((ASIN((refractive1/refractive2)*&
        (SIN(test_angle1*2*(3.14159265)/360))))*360/(2*3.14159265))
      END IF
WRITE(*,*) 'refract'
      refrac_num2=refrac_num1+1
    END IF

    IF (beam_angle2>=360) THEN          !This reduces the angle of the out going beam to an equivulent value between 0 and 360
      beam_angle2=beam_angle2-360
    END IF    

    WRITE(*,*) test_angle1
    WRITE(*,*) beam_angle2
    
END SUBROUTINE


SUBROUTINE true_intersection(old_intersect, beam_angle1, true_intersect, refractive, slope)
USE light_mods
IMPLICIT NONE

!This subroutine determines the true intersection

REAL, INTENT(IN)::old_intersect(2), beam_angle1
REAL, INTENT(OUT)::true_intersect(2), refractive, slope
REAL::test_length(1000), length, beam_angle2
INTEGER::i, j, ierror, num_range
TYPE (cross)::test_intersect(1000)

    OPEN(UNIT=14, FILE='prisms/test_intersects.f90', IOSTAT=ierror)
    
    i=0
    DO          !READ in the intersection locations and attached data
      i=i+1
       READ(14,*,IOSTAT=ierror) test_intersect(i)%refrac_index, test_intersect(i)%side(1:4), test_intersect(i)%intersect(1:2)
      IF (ierror /= 0) EXIT      !test for end of file   
    num_range=i          
    END DO

    IF (beam_angle1<0) THEN
      beam_angle2=360+beam_angle1
    ELSE
      beam_angle2=beam_angle1
    END IF

length=1000000    !This section kicks out any intersections that are in the wrong direction relative to the old intersection
    DO j=1, num_range
      IF ((beam_angle2<180) .AND. (beam_angle2>0) .AND. (test_intersect(j)%intersect(2)<old_intersect(2))) THEN
        WRITE(*,*) '1'
        CYCLE
      ELSE IF ((beam_angle2>180) .AND. (beam_angle2<360) .AND. (test_intersect(j)%intersect(2)>old_intersect(2))) THEN
        WRITE(*,*) '2'
        CYCLE
      ELSE IF (((beam_angle2==0) .OR. (beam_angle2==360)) .AND. (test_intersect(j)%intersect(1)<old_intersect(1))) THEN
        WRITE(*,*) '3'
        CYCLE
      ELSE IF ((beam_angle2==180) .AND. (test_intersect(j)%intersect(1)>old_intersect(1))) THEN
        WRITE(*,*) '4'
        CYCLE
      END IF

      !This section finds the closest intersection
      IF (((old_intersect(1)-test_intersect(j)%intersect(1))==0.0) .AND.&
      ((old_intersect(2)-test_intersect(j)%intersect(2))==0.0)) THEN
        CYCLE
      END IF
      test_length(j)=SQRT(((old_intersect(1)-test_intersect(j)%intersect(1))**2)+&
        ((old_intersect(2)-test_intersect(j)%intersect(2))**2))
      IF ((test_length(j) <= length) .AND. (test_length(j) >0.0001)) THEN
        length=test_length(j)
        true_intersect(1:2)=test_intersect(j)%intersect(1:2)
        refractive=test_intersect(j)%refrac_index
        IF (test_intersect(j)%side(1)==0) THEN
          slope=100000000
        ELSE
          slope=test_intersect(j)%side(2)
        END IF
      END IF

    END DO
    CLOSE(14)
    WRITE(*,*) true_intersect(1:2)
    WRITE(*,*) ''

    
END SUBROUTINE



SUBROUTINE intersect(source, test_surface, test_vertex, all_intersect)
IMPLICIT NONE

REAL, INTENT(IN):: source(4), test_surface(4), test_vertex(2,2)
REAL(KIND=4), INTENT(OUT):: all_intersect(2)
REAL(KIND=8)::source_test(3), test_surface2(3), test1(3), test2(3), dummyvar
REAL::vert(2,2)

!     This subroutine acts as a row reduction program to find the intersection of the light beam
!   and the the sides of a prisms.  The subroutine first converts the arrays of 4 to arrays of 3
!   of the form x y b to assist in the rref.

    IF (source(1)==0) THEN
      source_test(1)=-1        !If the source is a vertical line x=b
      source_test(2)=0
      source_test(3)=(source(4))
    ELSE 
      source_test(1)=(source(2))*(source(3))   !All source lines not of the form x=b
      source_test(2)=(source(1))
      source_test(3)=(source(4))
    END IF
        
    IF (test_surface(1)==0) THEN
      test_surface2(1)=-1+(source_test(1))            !If the test_surface is a vertical line x=b
      test_surface2(2)=0+(source_test(2))
      test_surface2(3)=(test_surface(4)+source(4))
    ELSE
      test_surface2(1)=(test_surface(2))*(test_surface(3))   !All test_surfaces lines not of the form x=b
      test_surface2(2)=(test_surface(1))
      test_surface2(3)=(test_surface(4))
    END IF
        
    dummyvar=(source_test(1))                               !Below is just solving of system of equations
    IF (dummyvar == 0) THEN
      dummyvar=(source_test(1)+test_surface2(1))
      test1(1)=source_test(1)+test_surface2(1)
      test1(2)=source_test(2)+test_surface2(2)
      test1(3)=source_test(3)+test_surface2(3)

    END IF
    test1(1)=(source_test(1))/(dummyvar)
    test1(2)=(source_test(2))/(dummyvar)
    test1(3)=(source_test(3))/(dummyvar)

!jc  array operations would be really nice here

    dummyvar=test_surface2(2)
    IF (dummyvar==0) THEN
      dummyvar=(source_test(2)+test_surface2(2))
      test2(1)=source_test(1)+test_surface2(1)
      test2(2)=source_test(2)+test_surface2(2)
      test2(3)=source_test(3)+test_surface2(3)
    END IF
    test2(1)=(test_surface2(1))/(dummyvar) 
    test2(2)=(test_surface2(2))/(dummyvar) 
    test2(3)=(test_surface2(3))/(dummyvar) 
    
    dummyvar=test1(2)
    test1(1)=(test1(1))-(test2(1))*(dummyvar)
    test1(2)=(test1(2))-(test2(2))*(dummyvar)
    test1(3)=(test1(3))-(test2(3))*(dummyvar)
    
    dummyvar=test1(1)
    test1(1)=(test1(1))/(dummyvar)
    test1(2)=(test1(2))/(dummyvar)
    test1(3)=(test1(3))/(dummyvar)
    
    dummyvar=test2(1)
    test2(1)=(test2(1))-(test1(1))*(dummyvar)
    test2(2)=(test2(2))-(test1(2))*(dummyvar)
    test2(3)=(test2(3))-(test1(3))*(dummyvar)
    
    test1(3)=(test1(3))*(-1)
    
    vert(1,1)=test_vertex(1,1)
    vert(1,2)=test_vertex(1,2)
    vert(2,1)=test_vertex(2,1)
    vert(2,2)=test_vertex(2,2)

    all_intersect(1)=test1(3)
    all_intersect(2)=test2(3)

    WRITE(*,*) test_vertex(1,1:2)
    WRITE(*,*) test_vertex(2,1:2)

    WRITE(*,*) all_intersect(1:2)

    WRITE(*,*) test1(1:3)
    WRITE(*,*) test2(1:3)

    
END SUBROUTINE


SUBROUTINE sides3(vertex, test_side)
IMPLICIT NONE

REAL, INTENT(IN)::vertex(3,2)
REAL, INTENT(OUT)::test_side(3,4)
REAL::test_vertex(4,2)
INTEGER::i

    !This subroutine creates three sided objects
    test_vertex(1,1:2)=vertex(1,1:2)
    test_vertex(2,1:2)=vertex(2,1:2)
    test_vertex(3,1:2)=vertex(3,1:2)    
    test_vertex(4,1:2)=vertex(1,1:2)

    DO i=1,3
      !Equation of test_side i  a*y=m*c*x+b
      !test_side(i,1)=a
      !test_side(i,2)=m
      !test_side(i,3)=c
      !test_side(i,4)=b
      IF ((test_vertex(i,2)-test_vertex((i+1),2)) == 0) THEN
        test_side(i,1)=1
        test_side(i,2)=0
        test_side(i,3)=0
        test_side(i,4)=(test_vertex(i,2))-((test_side(i,2))*(test_side(i,3))*(test_vertex(i,1)))
      ELSE IF ((test_vertex(i,1)-test_vertex((i+1),1)) == 0) THEN
        test_side(i,1)=0
        test_side(i,2)=1
        test_side(i,3)=1
        test_side(i,4)=test_vertex(i,1)
      ELSE
        test_side(i,1)=1
        test_side(i,2)=(test_vertex(i,2)-test_vertex((i+1),2))/(test_vertex(i,1)-test_vertex((i+1),1))
        test_side(i,3)=(test_vertex(i,1)-test_vertex((i+1),1))/(test_vertex(i,1)-test_vertex((i+1),1)) 
        test_side(i,4)=(test_vertex((i+1),2))-((test_side(i,2))*(test_side(i,3))*(test_vertex((i+1),1)))
      END IF
      WRITE(*,*) test_side(i,1:4)
    END DO
    
END SUBROUTINE


SUBROUTINE sides4(vertex, test_side)
IMPLICIT NONE

REAL, INTENT(IN)::vertex(4,2)
REAL, INTENT(OUT)::test_side(4,4)
REAL::test_vertex(5,2)
INTEGER::i

    !This subroutine creates four sided objects
    test_vertex(1,1:2)=vertex(1,1:2)
    test_vertex(2,1:2)=vertex(2,1:2)
    test_vertex(3,1:2)=vertex(3,1:2)    
    test_vertex(4,1:2)=vertex(4,1:2)
    test_vertex(5,1:2)=vertex(1,1:2)

    DO i=1,4
      !Equation of test_side i  a*y=m*c*x+b
      !test_side(i,1)=a
      !test_side(i,2)=m
      !test_side(i,3)=c
      !test_side(i,4)=b
      IF ((test_vertex(i,2)-test_vertex((i+1),2)) == 0) THEN
        test_side(i,1)=1
        test_side(i,2)=0
        test_side(i,3)=0
        test_side(i,4)=(test_vertex(i,2))-((test_side(i,2))*(test_side(i,3))*(test_vertex(i,1)))
      ELSE IF ((test_vertex(i,1)-test_vertex((i+1),1)) == 0) THEN
        test_side(i,1)=0
        test_side(i,2)=1
        test_side(i,3)=1
        test_side(i,4)=test_vertex(i,1)
      ELSE
        test_side(i,1)=1
        test_side(i,2)=(test_vertex(i,2)-test_vertex((i+1),2))/(test_vertex(i,1)-test_vertex((i+1),1))
        test_side(i,3)=(test_vertex(i,1)-test_vertex((i+1),1))/(test_vertex(i,1)-test_vertex((i+1),1)) 
        test_side(i,4)=(test_vertex((i+1),2))-((test_side(i,2))*(test_side(i,3))*(test_vertex((i+1),1)))
      END IF
    END DO

END SUBROUTINE


!jc  Very impressive program - it does a lot. My main issues with it are
!jc  that it could use more descriptive comments.  The size  and clarity 
!jc  of the code could have been reduced quite a bit by using array operations
!jc  and more subroutines.  Also, with all of the input that you take from
!jc  the screen, having more error checking and recovery from errors would
!jc  be good.  That said, you obviosuly put a lot of work into this.
