pro init_fc

common inp_block,fit_method,geometry,cp_file,sca_in,fc_init_file,sca_out,fc_out,fc_init,float_mode
common k_block,n_station,n_obs,n_sfo,n_stn,n_sfc,n_par
common l_block,l,stc,ssda,ftc,ftcr,ssdar,stcr,rfc,rsc,fbc,sbc,t,ti,s,si,ur,error,l_cnt
common i_block,stn_idx,fc_idx,mds,mrs,mdf,mrf,stca
common cp_block,point_id,numb_stn,cp_set,na,relaz,instr,n_st,n_cp
common ba_block,desired_stations
common ex_block,line_idx
common init_block,flng1,flat1
common planet_block,radius,lng0,lat0

; v2 6/24/4 include only those feature coordinates that have been floated

; v3 7/1/4 flat planet
;          Assume flat geometry, 0 feature altitude and use station coordinates
;          to initialize feature positions

; v4 spherical and flat planet
;      note that the main loop below counts every other floated feature coordinate which
;      assumes that z will be permanently non-floated and that fixed feature coordinates
;      are fixed in pairs, in other words, that all floated feature coordinates are also
;      floated in pairs.  YOU CANNOT USE THIS ROUTINE IF when fc X is floated and fc Y is not
;      also floated, or vice versa.  FEATURE COORDINATES MUST ALWAYS be FIXED or FLOATED
;      in PAIRS to use this routine.  As a token of this reality, floated feature coordinates
;      are indexed with m for both lng and lat, instead of m and m+1 for lng and lat
;         to derive initial feature coordinate guesses this routine substitutes
;         transposed elements of t for inverse t elements, i.e., t(0,1,vs(0)) instead of ti(1,0,vs(0)

; v4m 10/25/04 moves radius from common ba_block to common planet_block, formerly zero_block

; convert orientation angles from degrees to radians
  ssdar(0,*)=!dtor*ssda(0,*)
  ssdar(1,*)=!dtor*ssda(1,*)
  ssdar(2,*)=!dtor*ssda(2,*)

; compute transformation matrix T
  rr=ssdar(0,*)
  pr=ssdar(1,*)
  yr=ssdar(2,*)
  t(0,0,*)= cos(pr)*cos(yr)
  t(1,0,*)= cos(pr)*sin(yr)
  t(2,0,*)=-sin(pr)
  t(0,1,*)=-cos(rr)*sin(yr)+sin(rr)*sin(pr)*cos(yr)
  t(1,1,*)= cos(rr)*cos(yr)+sin(rr)*sin(pr)*sin(yr)
  t(2,1,*)= sin(rr)*cos(pr)
  t(0,2,*)= sin(rr)*sin(yr)+cos(rr)*sin(pr)*cos(yr)
  t(1,2,*)=-sin(rr)*cos(yr)+cos(rr)*sin(pr)*sin(yr)
  t(2,2,*)= cos(rr)*cos(pr)

  ftc(2,*)=0d0

  s_fc=size(fc_idx)
  if s_fc(0) ne 0 then begin
    for m=0,s_fc(1)-1,2 do begin
      ac= line_idx(mdf(m))
      ac1=line_idx(mdf(m+1))
      pid =point_id(ac)
      pid1=point_id(ac1)
      me=numb_stn(pid)-1

      xx_sum=0d0
      yy_sum=0d0
      sum_cnt=0d0

      for m1=0,me do begin
        vs=where(desired_stations eq cp_set(m1,pid),v_cnt)
        if v_cnt ne 0 then begin
          nar =!dtor*na(m1,pid)
          phar=!dtor*relaz(m1,pid)

          uu=t(0,0,vs(0))*sin(nar)*cos(phar)+t(0,1,vs(0))*sin(nar)*sin(phar)-t(0,2,vs(0))*cos(nar)
          vv=t(1,0,vs(0))*sin(nar)*cos(phar)+t(1,1,vs(0))*sin(nar)*sin(phar)-t(1,2,vs(0))*cos(nar)
          ww=t(2,0,vs(0))*sin(nar)*cos(phar)+t(2,1,vs(0))*sin(nar)*sin(phar)-t(2,2,vs(0))*cos(nar)

          if geometry eq 'f' then begin
            xx=stc(0,vs(0))-stc(2,vs(0))*uu/ww
            yy=stc(1,vs(0))-stc(2,vs(0))*vv/ww
            xx_sum=xx_sum+xx
            yy_sum=yy_sum+yy
          endif

          if geometry eq 's' then begin
            dellng=-!radeg*stc(2,vs(0))*uu/ww/radius
            dellat=-!radeg*stc(2,vs(0))*vv/ww/radius
            flat1(pid ,vs(0))= stc(1,vs(0))+dellat
            flng1(pid1,vs(0))= stc(0,vs(0))+dellng/cos(!dtor*flat1(pid,vs(0)))
            xx_sum=xx_sum+flng1(pid1,vs(0))
            yy_sum=yy_sum+flat1(pid ,vs(0))
          endif

          sum_cnt=sum_cnt+1d0

        endif
      endfor
      ftc(0,pid)=xx_sum/sum_cnt
      ftc(1,pid)=yy_sum/sum_cnt
;     ftc(point_id(line_idx(mrf(m  ))),pid )=xx_sum/sum_cnt
;     ftc(point_id(line_idx(mrf(m+1))),pid1)=yy_sum/sum_cnt
;     ftc(mrf(m  ),pid )=xx_sum/sum_cnt
;     ftc(mrf(m+1),pid1)=yy_sum/sum_cnt
    endfor
  endif
;stop
return
end