/[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.6 by mocchiut, Tue May 15 14:31:14 2012 UTC revision 1.10 by mocchiut, Fri Aug 29 07:06:41 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 607  c     COEF3 = COEFPATH(1:16)//COEF3 Line 607  c     COEF3 = COEFPATH(1:16)//COEF3
607  c     FILMOD(1) = COEF1  c     FILMOD(1) = COEF1
608  c     FILMOD(2) = COEF2  c     FILMOD(2) = COEF2
609  c     FILMOD(3) = COEF3  c     FILMOD(3) = COEF3
610        print *, "qui"  c      print *, "qui"
611        FILMOD(1) = P1(1:L1)        FILMOD(1) = P1(1:L1)
612        FILMOD(2) = P2(1:L2)        FILMOD(2) = P2(1:L2)
613        FILMOD(3) = P3(1:L3)        FILMOD(3) = P3(1:L3)
614        print *, "qua"  c      print *, "qua"
615  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'
616  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'
617  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'
618        WRITE(*,*) FILMOD(1)  c      WRITE(*,*) FILMOD(1)
619        WRITE(*,*) FILMOD(2)  c      WRITE(*,*) FILMOD(2)
620        WRITE(*,*) FILMOD(3)  c      WRITE(*,*) FILMOD(3)
621  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/
622        DATA   DTEMOD / 2005., 2010., 2015./        DATA   DTEMOD / 2005., 2010., 2015./
623  c      c    
# Line 633  C     ### numye = numye + 1 ; is number Line 633  C     ### numye = numye + 1 ; is number
633  C      C    
634  c     NUMYE=13  c     NUMYE=13
635        NUMYE=2        NUMYE=2
636        print *, "quo"  c      print *, "quo"
637                
638  C      C    
639  C     IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION  C     IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION
# Line 652  c     L = (IYEA - 1945)/5 + 1 Line 652  c     L = (IYEA - 1945)/5 + 1
652        FIL1 = FILMOD(L)          FIL1 = FILMOD(L)  
653        DTE2 = DTEMOD(L+1)        DTE2 = DTEMOD(L+1)
654        FIL2 = FILMOD(L+1)        FIL2 = FILMOD(L+1)
655        WRITE(*,*) FIL1  c      WRITE(*,*) FIL1
656        WRITE(*,*) FIL2  c      WRITE(*,*) FIL2
657        print *, "que"  c      print *, "que"
658  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
659        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
660        IF (IER .NE. 0) STOP                                  IF (IER .NE. 0) STOP                          
661        print *, "quessss"  c      print *, "quessss"
662        CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
663        IF (IER .NE. 0) STOP                            IF (IER .NE. 0) STOP                    
664        print *, "quj"  c      print *, "quj"
665  C--   DETERMINE IGRF COEFFICIENTS FOR YEAR  C--   DETERMINE IGRF COEFFICIENTS FOR YEAR
666        IF (L .LE. NUMYE-1) THEN                                IF (L .LE. NUMYE-1) THEN                        
667           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,
# Line 670  C--   DETERMINE IGRF COEFFICIENTS FOR YE Line 670  C--   DETERMINE IGRF COEFFICIENTS FOR YE
670           CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,               CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,    
671       1        GH2, NMAX, GHA)                                           1        GH2, NMAX, GHA)                                    
672        ENDIF        ENDIF
673        print *, "quw"  c      print *, "quw"
674  C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G  C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
675        F0=0.D0        F0=0.D0
676        DO 1234 J=1,3        DO 1234 J=1,3
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
# Line 744  C     ---------------------------------- Line 744  C     ----------------------------------
744        WRITE(FOUT,667) FSPEC        WRITE(FOUT,667) FSPEC
745  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)
746   667  FORMAT(A258)   667  FORMAT(A258)
747        print *," gui"  c      print *," gui"
748        OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)            OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)    
749        print *," gua"  c      print *," gua"
750        READ (IU, *, IOSTAT=IER, ERR=999)        READ (IU, *, IOSTAT=IER, ERR=999)
751        print *," gue"  c      print *," gue"
752        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD
753        print *," guo"  c      print *," guo"
754  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
755  C       Read the coefficient file, arranged as follows:                C       Read the coefficient file, arranged as follows:              
756  C                                                                                C                                                                              
# Line 787  C -------------------------------------- Line 787  C --------------------------------------
787                  ENDIF                                                                  ENDIF                                                
788  2233        CONTINUE                                                      2233        CONTINUE                                                    
789  2211    CONTINUE                                                          2211    CONTINUE                                                        
790        print *," guj"  c      print *," guj"
791                                                                                                                                                                    
792  999     CLOSE (IU)                                                    999     CLOSE (IU)                                                  
793        print *," guw IER",IER  c      print *," guw IER",IER
794        if ( IER .eq. -1 ) IER = 0 ! gfortran 4.1.2 bug workaround... hoping not to create problems with other versions        if ( IER .eq. -1 ) IER = 0 ! gfortran 4.1.2 bug workaround... hoping not to create problems with other versions
795                                                                                                                                                                    
796          RETURN                                                                RETURN                                                      
# Line 947  C ERA, EREQU and ERPOL as recommended by Line 947  C ERA, EREQU and ERPOL as recommended by
947  C ASTRONOMICAL UNION .  C ASTRONOMICAL UNION .
948  C-----------------------------------------------------------------  C-----------------------------------------------------------------
949          INTEGER TL1,TL2,TL3          INTEGER TL1,TL2,TL3
950          CHARACTER (len=258) TP1,TP2,TP3          CHARACTER (len=*) :: TP1,TP2,TP3
951          INTEGER L1,L2,L3          INTEGER L1,L2,L3
952          CHARACTER *258 P1,P2,P3          CHARACTER *258 P1,P2,P3
953          COMMON/PPATH/ L1,L2,L3,P1, P2, P3          COMMON/PPATH/ L1,L2,L3,P1, P2, P3
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.23