PROGRAM DBCOMP                                                    00001020
C
C     dbcomp_mars.f 
C
C      Calculates the vector (Br,Btheta,Bphi) fields
C      from an equivalent source magnetization distribution
C                                                                               
      parameter (ntot=11562,ipro=30000)
      Real slar(ipro),slor(ipro),salr(ipro)
      Real brr(ipro),btr(ipro),bpr(ipro)
C
      REAL P(ntot),CC(ntot),CS(ntot),SS(ntot),bfield(ntot)                     
      REAL FLAT(ntot),FLON(ntot),CTH(ntot),STH(ntot)   
C                                                                               
      PI = 3.1415927                                                            
      arc = 0.0174533
      r = 3393.5
      cutoff= 1500.
      t3 = 3.0
      t2 = 2.0
      scaling = 100.0
      itot = 0
      NP= 11550                                                                 
C
C     COMPUTE VOLUME FROM DIPOLE SPACING
C
      DSP = 1.89                                                                
      VOL= (59.2**2)*(DSP**2)*40.0  
C
C     Read in a profile 
C
      Open(20,
     &  file='/home/purucker/mgs_cd_read/night_passes/99073',
     &  status='old',form='formatted')
C
      i = 0
      i1 = 0
 702  Read(20,*,end=88) qt,sla,slo,sal,bx,by,bz
      i = i+1
C
C     Select a profile across Cimmeria
C
      If (i.le.18974.or.i.gt.22027) go to 702
C
      i1 = i1+1
      slar(i1)=sla
      slor(i1)=slo
      salr(i1)=sal
      ct=cos((90.-sla)*arc)
      st=sin((90.-sla)*arc)
      cp=cos(slo*arc)
      sp=sin(slo*arc)
      brr(i1)=bz*ct+(bx*cp+by*sp)*st
      btr(i1)=-bz*st+(bx*cp+by*sp)*ct
      bpr(i1)=by*cp-bx*sp
      Go to 702
  88  iproftot=i1
C                                                                               
C     global equivalent source solution based on night-time Br mapping
C     data, day and night Br aerobraking data, and night Br science
C     phasing orbit data 
C
      Open(15,
     & file='global.xyz',status='old',form='formatted')
C
      Do 18 i=1,np
      READ(15,92) FLON(I),FLAT(I),P(I),ss(i),cc(i),cs(i)
 92   Format(2F10.2,F15.4,3F10.3)
 18   CONTINUE                                                                  
C
C     Output file definition
C
      Open(40,file='cimmeria_mapping_profile.99073',
     &    status='new',form='formatted')
C 
      Do 510 m=1,iproftot
      xyz1=slar(m)
      xyz2=slor(m)
      xyz3=salr(m)
      T=R+XYZ3                                     
C                                                   
      TCR=(90.-XYZ1)*ARC                             
      CTC=COS(TCR)                                      
      STC=SIN(TCR)                                       
C                                                                               
      FR = 0.0                                                          
      FT = 0.0                                                          
      FP = 0.0                                                          
C                                                                               
C     SUM THE CONTRIBUTIONS FROM THE INDIVIDUAL DIPOLES                         
C                                                                               
      DO 57 I=1,NP                                         
        CTD=sin(flat(I)*arc)
        STD=cos(flat(i)*arc)
        DPH=(XYZ2-FLON(I))*ARC
        CDP=COS(DPH)
        SDP=SIN(DPH)
        DD=SS(I)
        BD=CC(I)
        CD=CS(I)
        A1=CTD*CTC+STD*STC*CDP
        B1=CTC*STD-STC*CTD*CDP
        C1=STC*SDP
        A2=-STC*CTD+CTC*STD*CDP
        B2=-STC*STD-CTC*CTD*CDP
        C2=CTC*SDP
        A3=-STD*SDP
        B3=CTD*SDP
        C3=CDP
        E = R*R + T*T - 2.0*R*T*A1
        C = SQRT(E)
        If (c.gt.1500.0) go to 57
        C = C * E
        C= (scaling*vol)/c
        F1= T*A1-R
        F2= -T*B1
        F3= T*C1
        D1=T-R*A1
        D2=-R*A2
        D3=-R*A3
        D1=3.*D1
        D2=3.*D2
        D3=3.*D3
        F1=F1/E
        F2=F2/E
        F3=F3/E
        D11=C*(D1*F1-A1)
        D12=C*(D1*F2+B1)
        D13=C*(D1*F3-C1)
        D21=C*(D2*F1-A2)
        D22=C*(D2*F2+B2)
        D23=C*(D2*F3-C2)
        D31=C*(D3*F1-A3)
        D32=C*(D3*F2+B3)
        D33=C*(D3*F3-C3)
C
        DFDP1=DD*D11+BD*D12+CD*D13
        DFDP2=DD*D21+BD*D22+CD*D23
        DFDP3=DD*D31+BD*D32+CD*D33
C
        FR=P(I)*DFDP1 + FR
        FT=P(I)*DFDP2 + FT
        FP=P(I)*DFDP3 + FP
   57   CONTINUE 
      Write(40,883) xyz1,xyz2,xyz3,fr,brr(m),ft,btr(m),fp,bpr(m)  
 883  Format(9F9.2)
 510  continue
C                                                                               
      End