! 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.1all_intersect(1))).OR. & ((vertex(j,1)+0.1>all_intersect(1)).AND.(vertex(j+1,1)-0.1all_intersect(2)).AND.(vertex(j+1,2)-0.1all_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)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(*,*) '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.