/[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.11 by mocchiut, Fri Oct 31 14:02:53 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 196  C--------------------------------------- Line 196  C---------------------------------------
196        SAVE /FIDB0/        SAVE /FIDB0/
197        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
198        SAVE /GENER/        SAVE /GENER/
199          REAL FLS
200  C  C
201  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3
202  C-- STEP IS STEP SIZE FOR FIELD LINE TRACING  C-- STEP IS STEP SIZE FOR FIELD LINE TRACING
# Line 203  C-- STEQ IS STEP SIZE FOR INTEGRATION Line 204  C-- STEQ IS STEP SIZE FOR INTEGRATION
204  C  C
205        DATA RMIN,RMAX    /0.05,1.01/        DATA RMIN,RMAX    /0.05,1.01/
206        DATA STEP,STEQ    /0.20,0.03/        DATA STEP,STEQ    /0.20,0.03/
207          BEQU=1.E10        BEQU=1.E10
208          FLS = FL
209  C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES  C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES
210        RLAT=GLAT*UMR        RLAT=GLAT*UMR
211        CT=SIN(RLAT)                                                      CT=SIN(RLAT)                                              
# Line 304  C*****INNER LOOP (FOR QUADRATURE)       Line 306  C*****INNER LOOP (FOR QUADRATURE)      
306        HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                                    HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                            
307        ZQ=Z*Z        ZQ=Z*Z
308        R=HLI+SQRT(HLI*HLI+ZQ)        R=HLI+SQRT(HLI*HLI+ZQ)
309          IF(R.NE.R)THEN
310             FL = FLS
311             RETURN
312          ENDIF
313        IF(R.LE.RMIN)GOTO30                                      IF(R.LE.RMIN)GOTO30                              
314        RQ=R*R        RQ=R*R
315        FF=SQRT(1.+3.*ZQ/RQ)                                      FF=SQRT(1.+3.*ZQ/RQ)                              
# Line 607  c     COEF3 = COEFPATH(1:16)//COEF3 Line 613  c     COEF3 = COEFPATH(1:16)//COEF3
613  c     FILMOD(1) = COEF1  c     FILMOD(1) = COEF1
614  c     FILMOD(2) = COEF2  c     FILMOD(2) = COEF2
615  c     FILMOD(3) = COEF3  c     FILMOD(3) = COEF3
616        print *, "qui"  c      print *, "qui"
617        FILMOD(1) = P1(1:L1)        FILMOD(1) = P1(1:L1)
618        FILMOD(2) = P2(1:L2)        FILMOD(2) = P2(1:L2)
619        FILMOD(3) = P3(1:L3)        FILMOD(3) = P3(1:L3)
620        print *, "qua"  c      print *, "qua"
621  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'
622  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'
623  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'
624        WRITE(*,*) FILMOD(1)  c      WRITE(*,*) FILMOD(1)
625        WRITE(*,*) FILMOD(2)  c      WRITE(*,*) FILMOD(2)
626        WRITE(*,*) FILMOD(3)  c      WRITE(*,*) FILMOD(3)
627  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/
628        DATA   DTEMOD / 2005., 2010., 2015./        DATA   DTEMOD / 2005., 2010., 2015./
629  c      c    
# Line 633  C     ### numye = numye + 1 ; is number Line 639  C     ### numye = numye + 1 ; is number
639  C      C    
640  c     NUMYE=13  c     NUMYE=13
641        NUMYE=2        NUMYE=2
642        print *, "quo"  c      print *, "quo"
643                
644  C      C    
645  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 658  c     L = (IYEA - 1945)/5 + 1
658        FIL1 = FILMOD(L)          FIL1 = FILMOD(L)  
659        DTE2 = DTEMOD(L+1)        DTE2 = DTEMOD(L+1)
660        FIL2 = FILMOD(L+1)        FIL2 = FILMOD(L+1)
661        WRITE(*,*) FIL1  c      WRITE(*,*) FIL1
662        WRITE(*,*) FIL2  c      WRITE(*,*) FIL2
663        print *, "que"  c      print *, "que"
664  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
665        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
666        IF (IER .NE. 0) STOP                                  IF (IER .NE. 0) STOP                          
667        print *, "quessss"  c      print *, "quessss"
668        CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
669        IF (IER .NE. 0) STOP                            IF (IER .NE. 0) STOP                    
670        print *, "quj"  c      print *, "quj"
671  C--   DETERMINE IGRF COEFFICIENTS FOR YEAR  C--   DETERMINE IGRF COEFFICIENTS FOR YEAR
672        IF (L .LE. NUMYE-1) THEN                                IF (L .LE. NUMYE-1) THEN                        
673           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,
# Line 670  C--   DETERMINE IGRF COEFFICIENTS FOR YE Line 676  C--   DETERMINE IGRF COEFFICIENTS FOR YE
676           CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,               CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,    
677       1        GH2, NMAX, GHA)                                           1        GH2, NMAX, GHA)                                    
678        ENDIF        ENDIF
679        print *, "quw"  c      print *, "quw"
680  C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G  C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
681        F0=0.D0        F0=0.D0
682        DO 1234 J=1,3        DO 1234 J=1,3
683           F = GHA(J) * 1.D-5           F = GHA(J) * 1.D-5
684           F0 = F0 + F * F           F0 = F0 + F * F
685   1234 CONTINUE   1234 CONTINUE
686        DIMO = DSQRT(F0)        DIMO = REAL(DSQRT(F0))
687                
688        GH1(1) =  0.0        GH1(1) =  0.0
689        I=2                  I=2          
# Line 693  c      print *, "quq" Line 699  c      print *, "quq"
699           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
700           F = F0 * 0.5D0                                               F = F0 * 0.5D0                                    
701           IF(IS.EQ.0) F = F * SQRT2           IF(IS.EQ.0) F = F * SQRT2
702           GH1(I) = GHA(I-1) * F0           GH1(I) = GHA(I-1) * REAL(F0)
703           I = I+1                                                   I = I+1                                        
704           DO 9 M=1,N                                               DO 9 M=1,N                                    
705              F = F * (X + M) / (X - M + 1.D0)                              F = F * (X + M) / (X - M + 1.D0)                
706              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))
707              GH1(I) = GHA(I-1) * F              GH1(I) = GHA(I-1) * REAL(F)
708              GH1(I+1) = GHA(I) * F              GH1(I+1) = GHA(I) * REAL(F)
709              I=I+2              I=I+2
710   9       CONTINUE                                             9       CONTINUE                                          
711        RETURN        RETURN
# Line 744  C     ---------------------------------- Line 750  C     ----------------------------------
750        WRITE(FOUT,667) FSPEC        WRITE(FOUT,667) FSPEC
751  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)
752   667  FORMAT(A258)   667  FORMAT(A258)
753        print *," gui"  c      print *," gui"
754        OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)            OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)    
755        print *," gua"  c      print *," gua"
756        READ (IU, *, IOSTAT=IER, ERR=999)        READ (IU, *, IOSTAT=IER, ERR=999)
757        print *," gue"  c      print *," gue"
758        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD
759        print *," guo"  c      print *," guo"
760  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
761  C       Read the coefficient file, arranged as follows:                C       Read the coefficient file, arranged as follows:              
762  C                                                                                C                                                                              
# Line 787  C -------------------------------------- Line 793  C --------------------------------------
793                  ENDIF                                                                  ENDIF                                                
794  2233        CONTINUE                                                      2233        CONTINUE                                                    
795  2211    CONTINUE                                                          2211    CONTINUE                                                        
796        print *," guj"  c      print *," guj"
797                                                                                                                                                                    
798  999     CLOSE (IU)                                                    999     CLOSE (IU)                                                  
799        print *," guw IER",IER  c      print *," guw IER",IER
800        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
801                                                                                                                                                                    
802          RETURN                                                                RETURN                                                      
# Line 947  C ERA, EREQU and ERPOL as recommended by Line 953  C ERA, EREQU and ERPOL as recommended by
953  C ASTRONOMICAL UNION .  C ASTRONOMICAL UNION .
954  C-----------------------------------------------------------------  C-----------------------------------------------------------------
955          INTEGER TL1,TL2,TL3          INTEGER TL1,TL2,TL3
956          CHARACTER (len=258) TP1,TP2,TP3          CHARACTER (len=*) :: TP1,TP2,TP3
957          INTEGER L1,L2,L3          INTEGER L1,L2,L3
958          CHARACTER *258 P1,P2,P3          CHARACTER *258 P1,P2,P3
959          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.11
  ViewVC Help
Powered by ViewVC 1.1.23