initialize_synthetic_camera Subroutine

public subroutine initialize_synthetic_camera(params, F)

Arguments

Type IntentOptional AttributesName
type(KORC_PARAMS), intent(in) :: params
type(FIELDS), intent(in) :: F

Contents


Source Code

  SUBROUTINE initialize_synthetic_camera(params,F)
    IMPLICIT NONE
    TYPE(KORC_PARAMS), INTENT(IN) :: params
    TYPE(FIELDS), INTENT(IN) :: F
    !	REAL(rp) :: aperture ! Aperture of the camera (diameter of lens)
    ! in meters
    REAL(rp) :: start_at
    !! in Seconds
    REAL(rp) :: Riw
    !! Radial position of inner wall
    INTEGER, DIMENSION(2) :: num_pixels
    !! Number of pixels (X,Y)
    REAL(rp), DIMENSION(2) :: sensor_size
    !! (horizontal,vertical)
    REAL(rp) :: focal_length
    REAL(rp), DIMENSION(2) :: position
    !! Position of camera (R,Z)
    REAL(rp) :: incline
    REAL(rp) :: lambda_min
    !! Minimum wavelength in cm
    REAL(rp) :: lambda_max
    !! Maximum wavelength in cm
    INTEGER :: Nlambda
    LOGICAL :: camera_on, integrated_opt, toroidal_sections, photon_count
    INTEGER :: ntor_sections
    REAL(rp) :: xmin, xmax, ymin, ymax, DX, DY
    INTEGER :: ii

    NAMELIST /SyntheticCamera/ camera_on,Riw,num_pixels,sensor_size,&
         focal_length,&
         position,incline,lambda_min,lambda_max,Nlambda,integrated_opt,&
         toroidal_sections,ntor_sections,start_at,photon_count

    open(unit=default_unit_open,file=TRIM(params%path_to_inputs),&
         status='OLD',form='formatted')
    read(default_unit_open,nml=SyntheticCamera)
    close(default_unit_open)

    !	write(*,nml=SyntheticCamera)

    cam%camera_on = camera_on

    if (cam%camera_on) then
       if (params%mpi_params%rank .EQ. 0) then
          write(6,'(/,"* * * * * * * * * * * * * * * * * *")')
          write(6,'("*  Initializing synthetic camera  *")')
       end if

       !	cam%aperture = aperture
       cam%start_at = start_at
       cam%Riw = Riw
       cam%num_pixels = num_pixels
       cam%sensor_size = sensor_size
       cam%pixel_area = PRODUCT(sensor_size/num_pixels)
       cam%focal_length = focal_length
       cam%position = position
       cam%incline = C_PI*incline/180.0_rp
       cam%horizontal_angle_view = ATAN2(0.5_rp*cam%sensor_size(1),&
            cam%focal_length)
       cam%vertical_angle_view = ATAN2(0.5_rp*cam%sensor_size(2),&
            cam%focal_length)
       cam%lambda_min = lambda_min ! In meters
       cam%lambda_max = lambda_max ! In meters
       cam%Nlambda = Nlambda
       cam%Dlambda = (cam%lambda_max - cam%lambda_min)/REAL(cam%Nlambda,rp)
       ALLOCATE(cam%lambda(cam%Nlambda))
       cam%photon_count = photon_count
       cam%integrated_opt = integrated_opt
       cam%toroidal_sections = toroidal_sections
       if (cam%toroidal_sections) then
          cam%ntor_sections = ntor_sections
       else
          cam%ntor_sections = 1_idef
       end if

       cam%r = (/COS(0.5_rp*C_PI - cam%incline),-SIN(0.5_rp*C_PI - &
            cam%incline),0.0_rp/)

       do ii=1_idef,cam%Nlambda
          cam%lambda(ii) = cam%lambda_min + REAL(ii-1_idef,rp)*cam%Dlambda
       end do

       ALLOCATE(cam%pixels_nodes_x(cam%num_pixels(1)))
       ALLOCATE(cam%pixels_nodes_y(cam%num_pixels(2)))
       ALLOCATE(cam%pixels_edges_x(cam%num_pixels(1) + 1))
       ALLOCATE(cam%pixels_edges_y(cam%num_pixels(2) + 1))

       xmin = -0.5_rp*cam%sensor_size(1)
       xmax = 0.5_rp*cam%sensor_size(1)
       DX = cam%sensor_size(1)/REAL(cam%num_pixels(1),rp)

       do ii=1_idef,cam%num_pixels(1)
          cam%pixels_nodes_x(ii) = xmin + 0.5_rp*DX + REAL(ii-1_idef,rp)*DX
       end do

       do ii=1_idef,cam%num_pixels(1)+1_idef
          cam%pixels_edges_x(ii) = xmin + REAL(ii-1_idef,rp)*DX
       end do

       ymin = -0.5_rp*cam%sensor_size(2)
       ymax = 0.5_rp*cam%sensor_size(2)
       DY = cam%sensor_size(2)/REAL(cam%num_pixels(2),rp)

       do ii=1_idef,cam%num_pixels(2)
          cam%pixels_nodes_y(ii) = ymin + 0.5_rp*DY + REAL(ii-1_idef,rp)*DY
       end do

       do ii=1_idef,cam%num_pixels(2)+1_idef
          cam%pixels_edges_y(ii) = ymin + REAL(ii-1_idef,rp)*DY
       end do

       ! Initialize ang variables
       ALLOCATE(ang%eta(cam%num_pixels(1)))
       ! angle between main line of sight and a given pixel along the x-axis
       ALLOCATE(ang%beta(cam%num_pixels(1))) 
       ALLOCATE(ang%psi(cam%num_pixels(2)+1_idef))
       ! angle between main line of sight and a given pixel along the y-axis

       do ii=1_idef,cam%num_pixels(1)
          ang%eta(ii) = ABS(ATAN2(cam%pixels_nodes_x(ii),cam%focal_length))
          if (cam%pixels_edges_x(ii) .LT. 0.0_rp) then
             ang%beta(ii) = 0.5_rp*C_PI - cam%incline - ang%eta(ii)
          else
             ang%beta(ii) = 0.5_rp*C_PI - cam%incline + ang%eta(ii)
          end if
       end do

       do ii=1_idef,cam%num_pixels(2)+1_idef
          ang%psi(ii) = ATAN2(cam%pixels_edges_y(ii),cam%focal_length)
       end do

       ang%threshold_angle = ATAN2(cam%Riw,-cam%position(1))
       ang%threshold_radius = SQRT(cam%Riw**2 + cam%position(1)**2)

       ALLOCATE(ang%ax(cam%num_pixels(1)))
       ! angle between main line of sight and a given pixel along the x-axis
       ALLOCATE(ang%ay(cam%num_pixels(2)))
       ! angle between main line of sight and a given pixel along the y-axis

       do ii=1_idef,cam%num_pixels(1)
          ang%ax(ii) = ATAN2(cam%pixels_nodes_x(ii),cam%focal_length)
       end do

       ang%min_ax = MINVAL(ang%ax)
       ang%max_ax = MAXVAL(ang%ax)

       do ii=1_idef,cam%num_pixels(2)
          ang%ay(ii) = ATAN2(cam%pixels_nodes_y(ii),cam%focal_length)
       end do

       ang%min_ay = MINVAL(ang%ay)
       ang%max_ay = MAXVAL(ang%ay)

       if (cam%incline.GT.0.5_rp*C_PI) then
          ang%ac = cam%incline - 0.5_rp*C_PI
       else
          ang%ac = 0.5_rp*C_PI - cam%incline
       end if


       ! Initialize poloidal plane parameters

       pplane%grid_dims = num_pixels
       ALLOCATE(pplane%nodes_R(pplane%grid_dims(1)))
       ALLOCATE(pplane%nodes_Z(pplane%grid_dims(2)))

       ! * * * * * * * ALL IN METERS * * * * * * * 

       IF (TRIM(params%field_model(1:10)) .EQ. 'ANALYTICAL') THEN
          pplane%Rmin = F%Ro - F%AB%a
          pplane%Rmax = F%Ro + F%AB%a
          pplane%Zmin = -F%AB%a
          pplane%Zmax = F%AB%a
       ELSE
          pplane%Rmin = MINVAL(F%X%R)
          pplane%Rmax = MAXVAL(F%X%R)
          pplane%Zmin = MINVAL(F%X%Z)
          pplane%Zmax = MAXVAL(F%X%Z)
       END IF

       pplane%DR = (pplane%Rmax - pplane%Rmin)/REAL(pplane%grid_dims(1),rp)
       pplane%DZ = (pplane%Zmax - pplane%Zmin)/REAL(pplane%grid_dims(2),rp)

       do ii=1_idef,pplane%grid_dims(1)
          pplane%nodes_R(ii) = pplane%Rmin + &
               0.5_rp*pplane%DR + REAL(ii-1_idef,rp)*pplane%DR
       end do

       do ii=1_idef,pplane%grid_dims(2)
          pplane%nodes_Z(ii) = pplane%Zmin + &
               0.5_rp*pplane%DZ + REAL(ii-1_idef,rp)*pplane%DZ
       end do
       ! * * * * * * * ALL IN METERS * * * * * * * 

       call save_synthetic_camera_params(params)

       if (params%mpi_params%rank .EQ. 0) then
          write(6,'("*     Synthetic camera ready!     *")')
          write(6,'("* * * * * * * * * * * * * * * * * *",/)')
       end if
    end if
  END SUBROUTINE initialize_synthetic_camera