cinitObs.f

Go to the documentation of this file.

!      cinitObs:  set up observation sites by reference to Zobsp.h
!        HeightList is obtained from DepthList.
!        However, if DepthList value is negative, HeightList has priority and
!        DepthList is obtained from HeightList.
!        Variables in Zobsv.h are calculated.
       subroutine cinitobs
       implicit none
!
#include  "Zglobalc.h"
#include  "Ztrackp.h"
#include  "Zcoord.h"
#include  "Zpos.h"
#include  "Zmagfield.h"
#include  "Zobs.h"
#include  "Zobsp.h"
#include  "Zobsv.h"
#include  "Zearth.h"
#include  "Zincidentp.h"
!
       integer i, icon
       logical error
       character*180  msg
       type(coord)::south 
       real*8  cvh2thick, cthick2h, mu
       external cblkobs
!
       noofsites = 0
       do i = 1, maxnoofsites
          if( depthlist(i) .eq. 0. ) goto 100    ! exit
          noofsites = noofsites + 1
          if( depthlist(i) .lt. 0. ) then
!                HeightList has priority
               depthlist(i) = cvh2thick(heightlist(i))
          else
               heightlist(i) = cthick2h(depthlist(i))
          endif
          call csetcoord('llh', latitofsite, 
     *    longitofsite, 
     *    heightlist(i),   obssites(i)%pos%xyz )
       enddo
 100   continue
       if(noofsites .eq. 0) then
          call cerrormsg('No of Observation sites = 0', 0)
       endif
!        set polar angle of the injection pointsite.
!        (used when ObsPlane==spherical)
       call csetcoord('llh', latitofsite, longitofsite, 
     *    heightofinj,   polarinjpos)
!        transform it to polar coord.
       call ctranscoord2('sph', polarinjpos, polarinjpos)

       if(endlevel .eq. 0) endlevel = noofsites   ! 96.09.18
       if(endlevel2 .eq. 0) endlevel2 = endlevel
!               set magnetic field for the deepest level. 
!                 Magfield is 'ned' here
       if(howgeomag .gt. 20 .and. howgeomag .le. 30 ) then
          magfieldned%x = magn
          magfieldned%y = mage
          magfieldned%z = magd
          magfieldned%sys = 'ned'
       else
!           if HowGeomag <=2, this is reset later. don't worry.
          call cgeomag(yearofgeomag, obssites(noofsites)%pos%xyz, 
     *              magfieldned, icon)
          if(icon .eq. 1) then
             write(msg, *) ' YearOfGeomag=',yearofgeomag, ' is ok ?'
             call cerrormsg(msg, 1)
          endif
       endif
!              to 'hva'
       call ctransmagto('hva',    obssites(noofsites)%pos%xyz, 
     * magfieldned, magfieldhva)
       call ctransmagto('xyz',    obssites(noofsites)%pos%xyz, 
     * magfieldned, magfieldxyz)
       if(abs(xaxisfromsouth) .gt. 360.d0) then
!           if it is > 360, use magnetic east to be the detector x-axis. 
!           angle form south to x-axis in counter closckwise.
!        Note:   'hva' system angle  + is clockwise.
!       XaxisFromSouth = 90.0 - MagfieldHVA.z
          xaxisfromsouth = 90.0 - magfieldhva%z
       endif

!          celestial coord. init.
       call kcelei(latitofsite, 
     *   longitofsite, dtgmt, 
     *   obssites(noofsites)%pos%height)
       call kadthi(xaxisfromsouth)

 
!          site coord is made  to be in the  xyz system.       
       do i = 1, noofsites
          call ctranscoord2('xyz', obssites(i)%pos%xyz,
     *        obssites(i)%pos%xyz )  ! to xyz
!                set ObsSites(i).pos.radiallen, ..depth, height
          call csetpos2(heightlist(i), depthlist(i), obssites(i)%pos)
          call cgetmoliereu(obssites(i)%pos%depth, 1.d0, mu)
          obssites(i)%mu = mu
       enddo
!          detector system coord. x is directed to mag. east.
!                                 z is directed to vertical.
!          get their direction cos. in 'xyz' system.
!
       coslatitude = cos(latitofsite*torad)
       sinlatitude = sin(latitofsite*torad)
       coslongitude = cos(longitofsite*torad)
       sinlongitude = sin(longitofsite*torad)
       detzaxis%r(1) = coslatitude * coslongitude
       detzaxis%r(2) = coslatitude * sinlongitude
       detzaxis%r(3) = sinlatitude
!          south direction is expressed as
       south%r(1) = sinlatitude *coslongitude
       south%r(2) = sinlatitude *sinlongitude
       south%r(3) = -coslatitude
!        in south-vertical system. X-axis is expressed
       detxaxis%r(1) = cos(xaxisfromsouth * torad)
       detxaxis%r(2) = sin(xaxisfromsouth * torad)
       detxaxis%r(3) = 0.
!           convert it to 'xyz' system       
       call ctransvectzx(1, detzaxis, south, detxaxis, detxaxis)
       txyz2det(1,:) = detxaxis%r(:)
       txyz2det(3,:) = detzaxis%r(:)
       call cvecprod(detzaxis%r(:), detxaxis%r(:), detyaxis%r(:))
       txyz2det(2,:) = detyaxis%r(:)
       tdet2xyz(:,:) = transpose(txyz2det(:,:))
!        check if in the ascending order of depth
       error = .false.
       do i = 2, noofsites
          if(depthlist(i-1) .ge. depthlist(i)) then
             write(msg, *) 'Depth is not in ascending order ',
     *       i, '-th one=',depthlist(i), 
     *       ' is <  previous one= ',  depthlist(i-1)
             call cerrormsg(msg, 1)
             error = .true.
          endif
       enddo
       if(error) then
          call cerrormsg('Obsev. depth/height not in order', 0)
       endif
       if(obsplane .eq. spherical .and.
     *       (za1ry .ne. 'is' .and. za1ry .ne. 'cos' ) ) then
          call cerrormsg(
     *      "ObsPlane==spherical but Za1ry !='is'/'cos'",0)
       endif
!       from v5.66
       if(borderheightl .eq. 0.) then
          borderheightl = obssites(noofsites)%pos%height -1.0
       endif
       obssites(0)%pos%height = borderheighth
       obssites(0)%pos%radiallen = borderheighth + eradius
       obssites(0)%pos%depth = cvh2thick(borderheighth)
       obssites(noofsites+1)%pos%height = borderheightl
       obssites(noofsites+1)%pos%radiallen = borderheightl + eradius
       obssites(noofsites+1)%pos%depth = cvh2thick(borderheightl)
       
!         +++++++++++++++++++++++ for AS ++++++++++++++++
       noofassites = 0
       do i = 1, maxnoofassites
          if( asdepthlist(i) .eq. 0. ) goto 200    ! exit
          noofassites = noofassites + 1
          if( asdepthlist(i) .lt. 0. ) then
!                HeightList has priority
               asdepthlist(i) = cvh2thick(asheightlist(i))
          else
               asheightlist(i) = cthick2h(asdepthlist(i))
          endif
          call csetcoord('llh', latitofsite, 
     *    longitofsite, 
     *    asheightlist(i),   asobssites(i)%pos%xyz )
       enddo
 200   continue
       
!          site coord is made  to be in the  xyz system.       
       do i = 1, noofassites
          call ctranscoord2('xyz', asobssites(i)%pos%xyz,
     *        asobssites(i)%pos%xyz )  ! to xyz
!                set ASObsSites(i).pos.radiallen, ..depth, height
          call csetpos2(asheightlist(i), asdepthlist(i), 
     *     asobssites(i)%pos)
!             Moliere Unit should be reset depending on
!             the 1ry angle. Here is a tentative one for vertical           
          call cgetmoliereu(asobssites(i)%pos%depth, 1.d0, mu)
          asobssites(i)%mu = mu
       enddo

       end