/[PAMELA software]/DarthVader/OrbitalInfo/src/igrf_sub.for
ViewVC logotype

Diff of /DarthVader/OrbitalInfo/src/igrf_sub.for

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.7 by mocchiut, Tue May 15 15:00:51 2012 UTC revision 1.9 by mocchiut, Fri Jan 17 09:54:31 2014 UTC
# Line 1  Line 1 
1          subroutine igrf_sub(xlat,xlong,year,height,  c        subroutine igrf_sub(xlat,xlong,year,height,
2       &          xl,icode,dip,dec)  c     &          xl,icode,dip,dec)
3  c----------------------------------------------------------------  c----------------------------------------------------------------
4  c   INPUT:  c   INPUT:
5  c     xlat     geodatic latitude in degrees  c     xlat     geodatic latitude in degrees
# Line 14  c          =3  an approximation is used Line 14  c          =3  an approximation is used
14  c     dip     geomagnetic inclination in degrees  c     dip     geomagnetic inclination in degrees
15  c     dec     geomagnetic declination in degress  c     dec     geomagnetic declination in degress
16  c----------------------------------------------------------------  c----------------------------------------------------------------
17    c
18        REAL              LATI,LONGI  c      REAL              LATI,LONGI
19        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD  c      COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
20        SAVE /GENER/  c      SAVE /GENER/
21  C  C
22        CALL INITIZE  c      CALL INITIZE
23        ibbb=0  c      ibbb=0
24        ALOG2=ALOG(2.)  c      ALOG2=ALOG(2.)
25        ISTART=1  c      ISTART=1
26          lati=xlat  c        lati=xlat
27          longi=xlong  c        longi=xlong
28  c  c
29  C----------------CALCULATE PROFILES-----------------------------------  C----------------CALCULATE PROFILES-----------------------------------
30  c  c
31          CALL FELDCOF(YEAR,DIMO)  c        CALL FELDCOF(YEAR,DIMO)
32          CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)  c        CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)
33          CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)  c        CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)
34          DIP=ASIN(BDOWN/BABS)/UMR  c        DIP=ASIN(BDOWN/BABS)/UMR
35        DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR  c      DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR
36        RETURN  c      RETURN
37        END  c      END
38  c  c
39  c  c
40  C SHELLIG.FOR, Version 2.0, January 1992  C SHELLIG.FOR, Version 2.0, January 1992
# Line 677  C--   DETERMINE MAGNETIC DIPOL MOMENT AN Line 677  C--   DETERMINE MAGNETIC DIPOL MOMENT AN
677           F = GHA(J) * 1.D-5           F = GHA(J) * 1.D-5
678           F0 = F0 + F * F           F0 = F0 + F * F
679   1234 CONTINUE   1234 CONTINUE
680        DIMO = DSQRT(F0)        DIMO = REAL(DSQRT(F0))
681                
682        GH1(1) =  0.0        GH1(1) =  0.0
683        I=2                  I=2          
# Line 693  c      print *, "quq" Line 693  c      print *, "quq"
693           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
694           F = F0 * 0.5D0                                               F = F0 * 0.5D0                                    
695           IF(IS.EQ.0) F = F * SQRT2           IF(IS.EQ.0) F = F * SQRT2
696           GH1(I) = GHA(I-1) * F0           GH1(I) = GHA(I-1) * REAL(F0)
697           I = I+1                                                   I = I+1                                        
698           DO 9 M=1,N                                               DO 9 M=1,N                                    
699              F = F * (X + M) / (X - M + 1.D0)                              F = F * (X + M) / (X - M + 1.D0)                
700              IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M))              IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M))
701              GH1(I) = GHA(I-1) * F              GH1(I) = GHA(I-1) * REAL(F)
702              GH1(I+1) = GHA(I) * F              GH1(I+1) = GHA(I) * REAL(F)
703              I=I+2              I=I+2
704   9       CONTINUE                                             9       CONTINUE                                          
705        RETURN        RETURN
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.23