C                  INTERNATIONAL AVS CENTER
C        (This disclaimer must remain at the top of all files)
C
C WARRANTY DISCLAIMER
C
C This module and the files associated with it are distributed free of charge.
C It is placed in the public domain and permission is granted for anyone to use,
C duplicate, modify, and redistribute it unless otherwise noted.  Some modules
C may be copyrighted.  You agree to abide by the conditions also included in
C the AVS Licensing Agreement, version 1.0, located in the main module
C directory located at the International AVS Center ftp site and to include
C the AVS Licensing Agreement when you distribute any files downloaded from
C that site.
C
C The International AVS Center, MCNC, the AVS Consortium and the individual
C submitting the module and files associated with said module provide absolutely
C NO WARRANTY OF ANY KIND with respect to this software.  The entire risk as to
C the quality and performance of this software is with the user.  IN NO EVENT
C WILL The International AVS Center, MCNC, the AVS Consortium and the individual
C submitting the module and files associated with said module BE LIABLE TO
C ANYONE FOR ANY DAMAGES ARISING FROM THE USE OF THIS SOFTWARE, INCLUDING,
C WITHOUT LIMITATION, DAMAGES RESULTING FROM LOST DATA OR LOST PROFITS, OR ANY
C SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES.
C
C This AVS module and associated files are public domain software unless
C otherwise noted.  Permission is hereby granted to do whatever you like with
C it, subject to the conditions that may exist in copyrighted materials. Should
C you wish to make a contribution toward the improvement, modification, or
C general performance of this module, please send us your comments:  why you
C liked or disliked it, how you use it, and most important, how it helps your
C work. We will receive your comments at avs@ncsc.org.
C
C Please send AVS module bug reports to avs@ncsc.org.
C
C This program simulates the reflection and refraction
C of electromagnetic waves from an interface.
C
C Written by David A Faux of the Physics Department,
C University of Surrey, England
C Used for teaching on Titan 1500.
C
C Revsions:
C   12 March 92  xfered via email to AVS UK.  - ianc
C   14 March 92  debugged  edit_properties call - ianc
C                moved filename to data block include
C
	INCLUDE '/usr/avs/include/avs.inc'
	INCLUDE '/usr/avs/include/geom.inc'
	INCLUDE 'waves.inc'
      external waves_desc
C   IAC CODE CHANGE :       include 'waves.dat'
	INCLUDE 'waves.dat'
C
C We initialize this module by passing a pointer to our description
C function.  The description function is defined below.  It defines our
C inputs, outputs and parameters.
C
      call AVScorout_init(waves_desc)
C
      loop_count = 0
      animate    = 0
C
      do while(.true.)
C
      input_result = AVScorout_input(displ_e_i, displ_h_i,
     &  displ_e_r, displ_h_r, displ_e_t, displ_h_t, animate,
     &  point(1), point(2), point(3), point(4), point(5),
     &  point(6), point(7), point(8), point(9), point(10),
     &  point(11), point(12),point(13), point(14),point(15),
     &  point(16), point(17),point(18), point(19),point(20),
     &  point(21) )
      call assign_vals
C
      if(animate.eq.0)then
      call update
      call AVScorout_wait()
      else
      call timestep
c*****
      call AVScorout_exec()
c*****
      endif
C
      loop_count = loop_count + 1
C
      enddo
C
      stop
      end
C
C This is the description routine called by AVScorout_init
C
      integer function waves_desc
	INCLUDE 'waves.inc'
	INCLUDE '/usr/avs/include/avs.inc'
	INCLUDE '/usr/avs/include/geom.inc'
C
C Set the module name and type
C
      call AVSset_module_name('waves demo', 'data')
C
C Create an output port for the result
C
      oport = AVScreate_output_port('Output Geometry', 'geom')
C
      parm = AVSadd_parameter('display Ei', 'boolean', 1, 0, 1)
      call AVSadd_parameter_prop(parm, 'width', 'integer', 2)
      call AVSadd_parameter_prop(parm, 'columns', 'integer', 2)
      parm = AVSadd_parameter('display Hi', 'boolean', 0, 0, 1)
      call AVSadd_parameter_prop(parm, 'width', 'integer', 2)
      call AVSadd_parameter_prop(parm, 'columns', 'integer', 2)
      parm = AVSadd_parameter('display Er', 'boolean', 0, 0, 1)
      call AVSadd_parameter_prop(parm, 'width', 'integer', 2)
      call AVSadd_parameter_prop(parm, 'columns', 'integer', 2)
      parm = AVSadd_parameter('display Hr', 'boolean', 0, 0, 1)
      call AVSadd_parameter_prop(parm, 'width', 'integer', 2)
      call AVSadd_parameter_prop(parm, 'columns', 'integer', 2)
      parm = AVSadd_parameter('display Et', 'boolean', 0, 0, 1)
      call AVSadd_parameter_prop(parm, 'width', 'integer', 2)
      call AVSadd_parameter_prop(parm, 'columns', 'integer', 2)
      parm = AVSadd_parameter('display Ht', 'boolean', 0, 0, 1)
C
      parm = AVSadd_parameter('Animate', 'boolean', 0, 0, 1)
      call AVSadd_parameter_prop(parm, 'width', 'integer', 2)
      call AVSadd_parameter_prop(parm, 'columns', 'integer', 2)
C
      rparm = AVSadd_parameter('E parallel', 'real', 1.0,
     &  0.0, 1.0)
      call typein_widget
      rparm = AVSadd_parameter('E perpendicular ', 'real', 0.0,
     &  0.0, 1.0)
      call typein_widget
      rparm = AVSadd_parameter('phase difference', 'real', 0.0,
     &  0.0, 90.0)
      call typein_widget
      rparm = AVSadd_parameter('amplify E      ', 'real', 0.2,
     &  0.0, 2.0)
      call typein_widget
      rparm = AVSadd_parameter('amplify H      ', 'real', 0.70,
     &  0.0, 2.0)
      call typein_widget
      rparm = AVSadd_parameter('Incident Angle ', 'real', 45.0,
     &  0.0, 90.0)
      call typein_widget
      rparm = AVSadd_parameter('refr. index - 1 ',
     &  'real', 1.0, 0.0, 20.0)
      call typein_widget
      rparm = AVSadd_parameter('refr. index - 2 ',
     &  'real', 1.5, 0.0, 20.0)
      call typein_widget
C
C The following are program outputs
C
      rparm = AVSadd_parameter('t_Ei_par', 'string',
     & '  ', '  ', '  ')
      call text_widget
      rparm = AVSadd_parameter('t_Er_par', 'string',
     & '  ', '  ', '  ')
      call text_widget
      rparm = AVSadd_parameter('t_Et_par', 'string',
     &'  ', '  ', '  ')
      call text_widget
      rparm = AVSadd_parameter('t_Ei_per', 'string',
     & '  ', '  ', '  ')
      call text_widget
      rparm = AVSadd_parameter('t_Er_per', 'string',
     &'  ', '  ', '  ')
      call text_widget
      rparm = AVSadd_parameter('t_Et_per', 'string',
     &   '  ', '  ', '  ')
      call text_widget
      rparm = AVSadd_parameter('t_p_par', 'string',
     &   '  ', '  ', '  ')
      call text_widget
      rparm = AVSadd_parameter('t_p_per', 'string',
     &   '  ', '  ', '  ')
      call text_widget
      rparm = AVSadd_parameter('t_ang_i', 'string',
     &'  ', '  ', '  ')
      call text_widget
      rparm = AVSadd_parameter('t_ang_r', 'string',
     &'  ', '  ', '  ')
      call text_widget
      rparm = AVSadd_parameter('t_ang_t', 'string',
     &'  ', '  ', '  ')
      call text_widget
      rparm = AVSadd_parameter('t_b_ang', 'string',
     &'  ', '  ', '  ')
      call text_widget
      rparm = AVSadd_parameter('t_c_ang', 'string',
     &'  ', '  ', '  ')
      call text_widget
C
      return
      end
C
C
      subroutine update
	INCLUDE '/usr/avs/include/avs.inc'
	INCLUDE '/usr/avs/include/geom.inc'
	INCLUDE 'waves.inc'
      integer i
C
      ri_rat = ri_1/ri_2
C
      call calc_angles
      call vertices
      if(loop_count.eq.0)call initial
C
      call filter
      call disp_verts
C
      call results
      call rescreen
C
      return
      end
C
C
      subroutine results
	INCLUDE '/usr/avs/include/avs.inc'
	INCLUDE '/usr/avs/include/geom.inc'
	INCLUDE 'waves.inc'
      write(cbuf,'(A,f6.4,A)')'Ei-par   ',amp_e_par,' V/m'
      call AVSmodify_parameter('t_Ei_par',AVS_VALUE,cbuf,0,0)
      write(cbuf,'(A,f6.4,A,f6.4,A)')'Er-par   ',f_r_par,
     & ' V/m  (R= ',R_par,')'
      call AVSmodify_parameter('t_Er_par',AVS_VALUE,cbuf,0,0)
      write(cbuf,'(A,f6.4,A,f6.4,A)')'Et-par   ',f_t_par,
     & ' V/m  (T= ',T_par,')'
      call AVSmodify_parameter('t_Et_par',AVS_VALUE,cbuf,0,0)
      write(cbuf,'(A,f6.4,A)')'Ei-per   ',amp_e_per,' V/m'
      call AVSmodify_parameter('t_Ei_per',AVS_VALUE,cbuf,0,0)
      write(cbuf,'(A,f6.4,A,f6.4,A)')'Er-per   ',f_r_per,
     & ' V/m  (R= ',R_per,')'
      call AVSmodify_parameter('t_Er_per',AVS_VALUE,cbuf,0,0)
      write(cbuf,'(A,f6.4,A,f6.4,A)')'Et-per   ',f_t_per,
     & ' V/m  (T= ',T_per,')'
      call AVSmodify_parameter('t_Et_per',AVS_VALUE,cbuf,0,0)
C
      write(cbuf,'(A,f7.2,A,f7.2)')'Phase-par (R) ',
     &(rad_to_deg*phi_r_par),' (T) ',(rad_to_deg*phi_t_par)
      call AVSmodify_parameter('t_p_par',AVS_VALUE,cbuf,0,0)
      write(cbuf,'(A,f7.2,A,f7.2)')'Phase-per (R) ',
     &(rad_to_deg*phi_r_per),' (T) ',(rad_to_deg*phi_t_per)
      call AVSmodify_parameter('t_p_per',AVS_VALUE,cbuf,0,0)
C
      write(cbuf,'(A,f8.4)')' Angle of i ',theta_i_deg
      call AVSmodify_parameter('t_ang_i',AVS_VALUE,cbuf,0,0)
      write(cbuf,'(A,f8.4)')' Angle of r ',theta_r_deg
      call AVSmodify_parameter('t_ang_r',AVS_VALUE,cbuf,0,0)
      write(cbuf,'(A,f8.4)')' Angle of t ',theta_t_deg
      call AVSmodify_parameter('t_ang_t',AVS_VALUE,cbuf,0,0)
      write(cbuf,'(A,f8.4)')' Brewster Angle ',b_angle
      call AVSmodify_parameter('t_b_ang',AVS_VALUE,cbuf,0,0)
      if(icrit.ge.1)then
      write(cbuf,'(A,f8.4)')' Critical Angle ',c_angle
      call AVSmodify_parameter('t_c_ang',AVS_VALUE,cbuf,0,0)
      else
      write(cbuf,'(A)')' Crit Angle    Not applicable '
      call AVSmodify_parameter('t_c_ang',AVS_VALUE,cbuf,0,0)
      endif
      return
      end
C
C
      subroutine vertices
C   IAC CODE CHANGE :       include 'waves.inc'
	INCLUDE 'waves.inc'
      real uxk, uyk, xh, yh, zh, distance, k_dot_r, dex, dey, dez
C
C first the lines for the e-wave, incident ray
C
      call fresnel
      call params
C
      do i=1,no_incr
      distance = real(i-1)*incr_1
      vert_e_i(1,1,i) = -distance*uxki
      vert_e_i(2,1,i) = -distance*uyki
      vert_e_i(3,1,i) = 0.0
      k_dot_r = 2.0*pi*(uxki*vert_e_i(1,1,i) +
     &   uyki*vert_e_i(2,1,i))/lambda_ru_1
      vert_e_i(1,2,i) = vert_e_i(1,1,i) + xei*cos(k_dot_r)
      vert_e_i(2,2,i) = vert_e_i(2,1,i) + yei*cos(k_dot_r)
      vert_e_i(3,2,i) = vert_e_i(3,1,i) + zei*cos(k_dot_r + phase_rad)
C
C ... and now the H-wave .. incident wave
C
      vert_h_i(1,1,i) = vert_e_i(1,1,i)
      vert_h_i(2,1,i) = vert_e_i(2,1,i)
      vert_h_i(3,1,i) = vert_e_i(3,1,i)
      dex = vert_e_i(1,2,i) - vert_e_i(1,1,i)
      dey = vert_e_i(2,2,i) - vert_e_i(2,1,i)
      dez = vert_e_i(3,2,i) - vert_e_i(3,1,i)
      call cross_prod(uxki, uyki, 0.0, dex, dey, dez, xh, yh, zh)
      vert_h_i(1,2,i) = vert_h_i(1,1,i) + xh*amp_h
      vert_h_i(2,2,i) = vert_h_i(2,1,i) + yh*amp_h
      vert_h_i(3,2,i) = vert_h_i(3,1,i) + zh*amp_h
C
      enddo
C
C Reflected ray
C
      do i=1,no_incr
      vert_e_r(1,1,i) = -vert_e_i(1,1,i)
      vert_e_r(2,1,i) = vert_e_i(2,1,i)
      vert_e_r(3,1,i) = 0.0
      k_dot_r = 2.0*pi*(uxkr*vert_e_r(1,1,i) +
     &   uykr*vert_e_r(2,1,i))/lambda_ru_1
      vert_e_r(1,2,i) = vert_e_r(1,1,i) + xer*cos(k_dot_r + phi_r_par)
      vert_e_r(2,2,i) = vert_e_r(2,1,i) + yer*cos(k_dot_r + phi_r_par)
      vert_e_r(3,2,i) = vert_e_r(3,1,i) + zer*cos(k_dot_r + phase_rad
     &   + phi_r_per)
C
C ... and now the H-wave .. reflected wave
C
      vert_h_r(1,1,i) = vert_e_r(1,1,i)
      vert_h_r(2,1,i) = vert_e_r(2,1,i)
      vert_h_r(3,1,i) = vert_e_r(3,1,i)
      dex = vert_e_r(1,2,i) - vert_e_r(1,1,i)
      dey = vert_e_r(2,2,i) - vert_e_r(2,1,i)
      dez = vert_e_r(3,2,i) - vert_e_r(3,1,i)
      call cross_prod(uxkr, uykr, 0.0, dex, dey, dez, xh, yh, zh)
      vert_h_r(1,2,i) = vert_h_r(1,1,i) + xh*amp_h
      vert_h_r(2,2,i) = vert_h_r(2,1,i) + yh*amp_h
      vert_h_r(3,2,i) = vert_h_r(3,1,i) + zh*amp_h
C
      enddo
C
C and the transmitted wave
C
      do i=1,no_incr
      distance = real(i-1)*incr_2
      vert_e_t(1,1,i) = distance*uxkt
      vert_e_t(2,1,i) = distance*uykt
      vert_e_t(3,1,i) = 0.0
      k_dot_r = 2.0*pi*(uxkt*vert_e_t(1,1,i) +
     &   uykt*vert_e_t(2,1,i))/lambda_ru_2
      vert_e_t(1,2,i) = vert_e_t(1,1,i) + xet*cos(k_dot_r + phi_t_par)
      vert_e_t(2,2,i) = vert_e_t(2,1,i) + yet*cos(k_dot_r + phi_t_par)
      vert_e_t(3,2,i) = vert_e_t(3,1,i) + zet*cos(k_dot_r + phase_rad
     &  + phi_t_per)
C
C ... and now the H-wave .. transmitted wave
C
      vert_h_t(1,1,i) = vert_e_t(1,1,i)
      vert_h_t(2,1,i) = vert_e_t(2,1,i)
      vert_h_t(3,1,i) = vert_e_t(3,1,i)
C
      if(icrit.le.1)then
      dex = vert_e_t(1,2,i) - vert_e_t(1,1,i)
      dey = vert_e_t(2,2,i) - vert_e_t(2,1,i)
      dez = vert_e_t(3,2,i) - vert_e_t(3,1,i)
      call cross_prod(uxkt, uykt, 0.0, dex, dey, dez, xh, yh, zh)
      vert_h_t(1,2,i) = vert_h_t(1,1,i) + xh*amp_h
      vert_h_t(2,2,i) = vert_h_t(2,1,i) + yh*amp_h
      vert_h_t(3,2,i) = vert_h_t(3,1,i) + zh*amp_h
C
      elseif(icrit.eq.2)then
C
      sin0t = ri_rat*sin(theta_i_rad)
      cos0t = sqrt(sin0t**2 - 1.0)
      vert_h_t(1,2,i) = vert_h_t(1,1,i) + zet*cos0t*
     &  cos(k_dot_r + phi_t_per - pi/2.0)
      vert_h_t(2,2,i) = vert_h_t(2,1,i) + zet*sin0t*
     &  cos(k_dot_r + phi_t_per)
      vert_h_t(3,2,i) = vert_h_t(3,1,i) - yet*sin0t*
     &  cos(k_dot_r + phi_t_par)  +  xet*cos0t*
     &  cos(k_dot_r + phi_t_per + pi/2.0)
C
      endif
C
      enddo
C
C Colors of waves
C
      call cols(ntab_i,col_e_i,col_e_i_0,col_e_i_m)
      call cols(ntab_i,col_h_i,col_h_i_0,col_h_i_m)
      call cols(ntab_r,col_e_r,col_e_r_0,col_e_r_m)
      call cols(ntab_r,col_h_r,col_h_r_0,col_h_r_m)
      call cols(ntab_t,col_e_t,col_e_t_0,col_e_t_m)
      call cols(ntab_t,col_h_t,col_h_t_0,col_h_t_m)
C
C Lines indicating wave trajectories
C
      do i=1,3
      vi(i,1) = vert_e_i(i,1,1)
      vi(i,2) = vert_e_i(i,1,no_incr)
      ci(i,1) = col_e_i_0(i)
      ci(i,2) = col_e_i_0(i)
      vr(i,1) = vert_e_r(i,1,1)
      vr(i,2) = vert_e_r(i,1,no_incr)
      cr(i,1) = col_e_r_0(i)
      cr(i,2) = col_e_r_0(i)
      vt(i,1) = vert_e_t(i,1,1)
      vt(i,2) = vert_e_t(i,1,no_incr)
      ct(i,1) = col_e_t_0(i)
      ct(i,2) = col_e_t_0(i)
      enddo
C
      return
      end
C
C
      subroutine cols(ntab,col,col0,colm)
C   IAC CODE CHANGE :       include 'waves.inc'
	INCLUDE 'waves.inc'
      real col(3,2,200), col0(3), colm(3)
      integer nval, i
C
      do i=1,no_incr
      if(mod(i,(no_lines/2)).eq.ntab)then
      col(1,1,i) = colm(1)
      col(2,1,i) = colm(2)
      col(3,1,i) = colm(3)
      col(1,2,i) = colm(1)
      col(2,2,i) = colm(2)
      col(3,2,i) = colm(3)
      else
      col(1,1,i) = col0(1)
      col(2,1,i) = col0(2)
      col(3,1,i) = col0(3)
      col(1,2,i) = col0(1)
      col(2,2,i) = col0(2)
      col(3,2,i) = col0(3)
      endif
      enddo
C
      return
      end
C
C
      subroutine calc_angles
C   IAC CODE CHANGE :       include 'waves.inc'
	INCLUDE 'waves.inc'
      theta_i_rad = deg_to_rad*theta_i_deg
      theta_r_rad = theta_i_rad
      theta_r_deg = rad_to_deg*theta_r_rad
      b_angle = rad_to_deg*atan(1./ri_rat)
      if(ri_rat.gt.1.0)then
      c_angle = rad_to_deg*asin(1./ri_rat)
        if(theta_i_deg.ge.c_angle)then
        theta_t_rad = pi/2.0
        icrit=2
        else
        theta_t_rad = asin(ri_rat*sin(theta_i_rad))
        icrit=1
        endif
      else
      theta_t_rad = asin(ri_rat*sin(theta_i_rad))
      icrit=0
      endif
      theta_t_deg = rad_to_deg*theta_t_rad
      phase_rad = deg_to_rad*phase_deg
      return
      end
C
C
      subroutine cross_prod(x1, y1, z1, x2, y2, z2, xa, ya, za)
C   IAC CODE CHANGE :       include 'waves.inc'
	INCLUDE 'waves.inc'
      real x1, y1, z1, x2, y2, z2, xa, ya, za, magk
      magk = sqrt(x1**2 + y1**2 + z1**2)
      xa = (y1*z2 - z1*y2)/2.0
      ya = (z1*x2 - x1*z2)/2.0
      za = (x1*y2 - y1*x2)/2.0
      return
      end
C
C
      real function value(a)
      real a
      value = a
      return
      end
C
C
      subroutine assign_vals
C   IAC CODE CHANGE :       include 'waves.inc'
	INCLUDE 'waves.inc'
      amp_e_par = value(%val(point(1)))
      amp_e_per = value(%val(point(2)))
      phase_deg = value(%val(point(3)))
      amp_e = value(%val(point(4)))
      amp_h = value(%val(point(5)))
      theta_i_deg = value(%val(point(6)))
      ri_1 = value(%val(point(7)))
      ri_2 = value(%val(point(8)))
      t_ei_par = value(%val(point(9)))
      t_er_par = value(%val(point(10)))
      t_et_par = value(%val(point(11)))
      t_ei_per = value(%val(point(12)))
      t_er_per = value(%val(point(13)))
      t_et_per = value(%val(point(14)))
      t_p_par = value(%val(point(15)))
      t_p_per = value(%val(point(16)))
      t_th_i_deg = value(%val(point(17)))
      t_th_r_deg = value(%val(point(18)))
      t_th_t_deg = value(%val(point(19)))
      t_b_ang = value(%val(point(20)))
      t_c_ang = value(%val(point(21)))
      return
      end
C
C
      subroutine fresnel
C   IAC CODE CHANGE :       include 'waves.inc'
	INCLUDE 'waves.inc'
      real u, v
      if(icrit.le.1)then
      u = cos(theta_i_rad)
      v = ri_rat*cos(theta_t_rad)
      f_r_par = (-u+v)/(u+v)
      f_t_par = 2.0*ri_rat*u/(u+v)
      if((-u+v).lt.0.0)then
        phi_r_par = pi
      else
        phi_r_par = 0.0
      endif
      phi_t_par = 0.0
      u = u*ri_rat
      v = v/ri_rat
      f_r_per = (u-v)/(u+v)
      f_t_per = 2.0*u/(u+v)
      if((u-v).lt.0.0)then
        phi_r_per = pi
      else
        phi_r_per = 0.0
      endif
      phi_t_per = 0.0
      R_par = f_r_par**2
      T_par = (cos(theta_t_rad)*f_t_par**2)
     &  /cos(theta_i_rad)/ri_rat
      R_per = f_r_per**2
      T_per = (cos(theta_t_rad)*f_t_per**2)
     &  /cos(theta_i_rad)/ri_rat
      f_r_par = abs(f_r_par*amp_e_par)
      f_r_per = abs(f_r_per*amp_e_per)
      f_t_par = abs(f_t_par*amp_e_par)
      f_t_per = abs(f_t_per*amp_e_per)
      if(amp_e_par.eq.0.0)then
        R_par = 0.0
        T_par = 0.0
      endif
      if(amp_e_per.eq.0.0)then
        R_per = 0.0
        T_per = 0.0
      endif
C
      else
C
      f_r_par = amp_e_par
      u = 2.0/sqrt(ri_rat**(-2) + 1)
      f_t_par = amp_e_par*u
      f_r_per = amp_e_per
      u = 2.0*cos(theta_i_rad)/sqrt(1-ri_rat**(-2))
      f_t_per = amp_e_per*u
      R_par = 1.0
      R_per = 1.0
      T_par = 0.0
      T_per = 0.0
      if(amp_e_par.eq.0.0)R_par = 0.0
      if(amp_e_per.eq.0.0)R_per = 0.0
      v = ri_rat*sin(theta_i_rad)
      phi_r_par = -2.0*atan(cos(theta_i_rad)/ri_rat**2/
     &  sqrt(v**2 - ri_rat**(-2)))
      phi_t_par = phi_r_par/2.0 - pi/2.
      phi_r_per = 2.0*atan(sqrt(v**2 - 1.0)/(ri_rat*cos(theta_i_rad)))
      phi_t_per = phi_r_per/2.0
      endif
      return
      end
C
C
      subroutine params
C   IAC CODE CHANGE :       include 'waves.inc'
	INCLUDE 'waves.inc'
      lambda_ru_1 = speed_light/freq/unit_length
      lambda_ru_2 = lambda_ru_1*ri_1/ri_2
      no_incr = no_lines*int(0.5/lambda_ru_1 + 0.5)
      incr_1 = lambda_ru_1/real(no_lines)
      incr_2 = lambda_ru_2/real(no_lines)
      uxki = -sin(theta_i_rad)
      uyki = cos(theta_i_rad)
      uxkr = uxki
      uykr = -uyki
      uxkt = -sin(theta_t_rad)
      uykt = cos(theta_t_rad)
      xei = -amp_e*amp_e_par*cos(theta_i_rad)
      yei = -amp_e*amp_e_par*sin(theta_i_rad)
      zei = amp_e*amp_e_per
      xer = xei*f_r_par
      yer = -yei*f_r_par
      zer = zei*f_r_per
      xet = -f_t_par*amp_e*amp_e_par*cos(theta_t_rad)
      yet = -f_t_par*amp_e*amp_e_par*sin(theta_t_rad)
      zet = zei*f_t_per
      return
      end
C
C
      subroutine initial
	INCLUDE '/usr/avs/include/avs.inc'
	INCLUDE '/usr/avs/include/geom.inc'
	INCLUDE 'waves.inc'
      integer i,j,k,nvs
      open (unit=10, file=filename)
      rewind 10
      do i=1,4
      intc2(1,i) = 1.0
      intc2(2,i) = 0.0
      intc2(3,i) = 0.0
      intc1(1,i) = 0.0
      intc1(2,i) = 0.0
      intc1(3,i) = 1.0
      enddo
C
      objn = geom_create_obj(GEOM_POLYTRI, GEOM_NULL)
      call geom_add_disjoint_line(objn,vn,cn,2,GEOM_COPY_DATA)
C
      objint1 = geom_create_obj(GEOM_POLYHEDRON, GEOM_NULL)
      objint2 = geom_create_obj(GEOM_POLYHEDRON, GEOM_NULL)
C
      do k=1,6
        read(10,*)nvs
        read(10,*)((intv(i,j),i=1,3), j=1,nvs)
        call geom_add_disjoint_polygon(objint1, intv, GEOM_NULL,
     &         intc1, nvs, GEOM_NOT_SHARED, 0)
        do j=1,nvs
            intv(2,j) = intv(2,j) + 0.5
        enddo
        call geom_add_disjoint_polygon(objint2, intv, GEOM_NULL,
     &         intc2, nvs, GEOM_NOT_SHARED, 0)
      enddo
C
      call geom_gen_normals(objint1,0)
      call geom_gen_normals(objint2,0)
      call geom_cvt_polyh_to_polytri(objint1,
     &  ior(ior(GEOM_SURFACE, GEOM_WIREFRAME),0))
      call geom_cvt_polyh_to_polytri(objint2,
     &  ior(ior(GEOM_SURFACE, GEOM_WIREFRAME),0))
C
      return
      end
C
C
      subroutine rescreen
	INCLUDE '/usr/avs/include/avs.inc'
	INCLUDE '/usr/avs/include/geom.inc'
	INCLUDE 'waves.inc'

      editlist = 0
      editlist = geom_init_edit_list(editlist)
      call geom_edit_geometry(editlist,'waves_demo',obj)
      call geom_edit_geometry(editlist,'waves_demo',objn)
c
c --- seems to be a bug in edit properties, so leave to the network file
c --- original
c      call geom_edit_geometry(editlist,'waves1',objint1)
c      call geom_edit_properties(editlist, 'waves1', -1.0, -1.0,
c     &  -1.0, -1.0, tint1, -1.0, -1.0, -1.0)
c      call geom_edit_geometry(editlist,'waves2',objint2)
c      call geom_edit_properties(editlist, 'waves2', -1.0, -1.0,
c     &  -1.0, -1.0, tint2, -1.0, -1.0, -1.0)
c --- new
      call geom_edit_geometry(editlist,'waves1',objint1)
      call geom_edit_properties(editlist, 'waves1', -1.0, -1.0,
     &  -1.0, -1.0, tint1, -1.0)
      call geom_edit_geometry(editlist,'waves2',objint2)
      call geom_edit_properties(editlist, 'waves2', -1.0, -1.0,
     &  -1.0, -1.0, tint2, -1.0)
c
      call AVScorout_output(editlist)
      call geom_destroy_edit_list(editlist)
      call geom_destroy_obj(obj)
      return
      end
C
C
      subroutine timestep
	INCLUDE '/usr/avs/include/avs.inc'
	INCLUDE '/usr/avs/include/geom.inc'
	INCLUDE 'waves.inc'
      istep = istep + 1
      if(istep.gt.(no_incr*ntsteps))istep = istep -
     &   (no_incr*ntsteps)
      call filter
      call disp_verts
      call rescreen
      return
      end
C
C
      subroutine disp_verts
	INCLUDE '/usr/avs/include/avs.inc'
	INCLUDE '/usr/avs/include/geom.inc'
	INCLUDE 'waves.inc'
C
      obj = geom_create_obj(GEOM_POLYTRI, GEOM_NULL)
      tint1 = (ri_1 - 1.0)/3.0
      if(tint1.gt.0.33)tint1=0.33
      tint2 = (ri_2 - 1.0)/3.0
      if(tint2.gt.0.33)tint2=0.33
C
      if(displ_e_i.eq.1)
     &  call geom_add_disjoint_line(obj,svert_e_i,col_e_i,
     &  (2*no_incr), GEOM_COPY_DATA)
      if(displ_e_r.eq.1)
     &  call geom_add_disjoint_line(obj,svert_e_r,col_e_r,
     &  (2*no_incr), GEOM_COPY_DATA)
      if(displ_e_t.eq.1)
     &  call geom_add_disjoint_line(obj,svert_e_t,col_e_t,
     &  (2*no_incr), GEOM_COPY_DATA)
      if(displ_h_i.eq.1)
     &  call geom_add_disjoint_line(obj,svert_h_i,col_h_i,
     &  (2*no_incr), GEOM_COPY_DATA)
      if(displ_h_r.eq.1)
     &  call geom_add_disjoint_line(obj,svert_h_r,col_h_r,
     &  (2*no_incr), GEOM_COPY_DATA)
      if(displ_h_t.eq.1)
     &  call geom_add_disjoint_line(obj,svert_h_t,col_h_t,
     &  (2*no_incr), GEOM_COPY_DATA)
C
      call geom_add_disjoint_line(obj,vi,ci,2,GEOM_COPY_DATA)
      call geom_add_disjoint_line(obj,vr,cr,2,GEOM_COPY_DATA)
      call geom_add_disjoint_line(obj,vt,ct,2,GEOM_COPY_DATA)
C
      return
      end
C
C
      subroutine filter
C   IAC CODE CHANGE :       include 'waves.inc'
	INCLUDE 'waves.inc'
      integer i,j,nadd
      real di(3), dr(3), dt(3)
      di(1) = real(istep)*incr_1*uxki/real(ntsteps)
      di(2) = real(istep)*incr_1*uyki/real(ntsteps)
      di(3) = 0.0
      dr(1) = real(istep)*incr_1*uxkr/real(ntsteps)
      dr(2) = real(istep)*incr_1*uykr/real(ntsteps)
      dr(3) = 0.0
      dt(1) = real(istep)*incr_2*uxkt/real(ntsteps)
      dt(2) = real(istep)*incr_2*uykt/real(ntsteps)
      dt(3) = 0.0
      do i=1,no_incr
        do j=1,3
        svert_e_i(j,1,i) = vert_e_i(j,1,i) + di(j)
        svert_e_i(j,2,i) = vert_e_i(j,2,i) + di(j)
        svert_h_i(j,1,i) = vert_h_i(j,1,i) + di(j)
        svert_h_i(j,2,i) = vert_h_i(j,2,i) + di(j)
        svert_e_r(j,1,i) = vert_e_r(j,1,i) + dr(j)
        svert_e_r(j,2,i) = vert_e_r(j,2,i) + dr(j)
        svert_h_r(j,1,i) = vert_h_r(j,1,i) + dr(j)
        svert_h_r(j,2,i) = vert_h_r(j,2,i) + dr(j)
        svert_e_t(j,1,i) = vert_e_t(j,1,i) + dt(j)
        svert_e_t(j,2,i) = vert_e_t(j,2,i) + dt(j)
        svert_h_t(j,1,i) = vert_h_t(j,1,i) + dt(j)
        svert_h_t(j,2,i) = vert_h_t(j,2,i) + dt(j)
        enddo
      enddo
      di(1) = no_incr*incr_1*uxki
      di(2) = no_incr*incr_1*uyki
      dr(1) = no_incr*incr_1*uxkr
      dr(2) = no_incr*incr_1*uykr
      dt(1) = no_incr*incr_2*uxkt
      dt(2) = no_incr*incr_2*uykt
      nadd = (istep+ntsteps-1)/ntsteps
      do i=1,nadd
        do j=1,2
        svert_e_i(j,1,i) = svert_e_i(j,1,i) - di(j)
        svert_e_i(j,2,i) = svert_e_i(j,2,i) - di(j)
        svert_h_i(j,1,i) = svert_h_i(j,1,i) - di(j)
        svert_h_i(j,2,i) = svert_h_i(j,2,i) - di(j)
        enddo
      enddo
      do i=no_incr,no_incr-nadd+2,-1
        do j=1,2
        svert_e_r(j,1,i) = svert_e_r(j,1,i) - dr(j)
        svert_e_r(j,2,i) = svert_e_r(j,2,i) - dr(j)
        svert_h_r(j,1,i) = svert_h_r(j,1,i) - dr(j)
        svert_h_r(j,2,i) = svert_h_r(j,2,i) - dr(j)
        svert_e_t(j,1,i) = svert_e_t(j,1,i) - dt(j)
        svert_e_t(j,2,i) = svert_e_t(j,2,i) - dt(j)
        svert_h_t(j,1,i) = svert_h_t(j,1,i) - dt(j)
        svert_h_t(j,2,i) = svert_h_t(j,2,i) - dt(j)
        enddo
      enddo
C
      return
      end
C
C
      subroutine text_widget
	INCLUDE '/usr/avs/include/avs.inc'
	INCLUDE '/usr/avs/include/geom.inc'
C   IAC CODE CHANGE :       include 'waves.inc'
	INCLUDE 'waves.inc'
      parm = AVSconnect_widget(rparm, 'text')
      call AVSadd_parameter_prop(rparm, 'width', 'integer', 4)
      call AVSadd_parameter_prop(rparm, 'columns', 'integer', 1)
      return
      end
C
C
      subroutine typein_widget
	INCLUDE '/usr/avs/include/avs.inc'
	INCLUDE '/usr/avs/include/geom.inc'
	INCLUDE 'waves.inc'
      parm = AVSconnect_widget(rparm, 'typein_real')
      call AVSadd_parameter_prop(rparm, 'width', 'integer', 4)
      call AVSadd_parameter_prop(rparm, 'columns', 'integer', 1)
      return
      end
C