text:
      SUBROUTINE LOCMIN(M,X,N,P,SSQ,ALF,EFS,ITRAP,ESCF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      COMMON /NLLSQI/ NCOUNT
      COMMON /KEYWRD/ CONTRL
      COMMON /NUMCAL/ NUMCAL
      DIMENSION X(*), P(*), EFS(*)
      DIMENSION CONST(MAXPAR), XSTOR(MAXPAR), GSTOR(MAXPAR)
      DIMENSION PHI(3),VT(3)
      INTEGER LEFT,RIGHT,CENTER
      CHARACTER*241 KEYWRD, CONTRL
      SAVE XMAXM, SCALE,  KEYWRD, EPS, DEBUG, TEE, YMAXST, XCRIT
      SAVE MXCNT2, IPRINT, CONST
      LOGICAL  DEBUG, LOWER
      DATA CONST/MAXPAR*1.D0/
      DATA ICALCN/0/
************************************************************************
*
*    LOCMIN IS CALLED BY NLLSQ ONLY. IT IS A LINE-SEARCH PROCEDURE FOR
*    LOCATING A MINIMUM IN THE FUNCTION SPACE OF COMPFG.  SEE NLLSQ
*    FOR MORE DETAILS
*
************************************************************************
      IF (ICALCN.NE.NUMCAL) THEN
         ICALCN=NUMCAL
         XMAXM=1.D9
         SCALE=1.D0
         KEYWRD=CONTRL
C
C THE ABOVE LINE IS TO TRY TO PREVENT OVERFLOW IN NLLSQ
C
         EPS=1.D-5
         DEBUG=(INDEX(KEYWRD,'LINMIN') .NE. 0)
         TEE=1.D-2
         YMAXST=0.005D0
         XCRIT=0.0002D0
         MXCNT2=30
         IPRINT=0
         IF(DEBUG)IPRINT=-1
      ENDIF
      XMAXM=1.D-11
      DO 10 I=1,N
   10 XMAXM=MAX(XMAXM,ABS(P(I)))
      XMINM=XMAXM*SCALE
      XMAXM=YMAXST/XMAXM/SCALE
      FIN = SSQ
      LOWER = .FALSE.
      T=ALF
      PHI(1) = SSQ
      VT(1) = 0.0D0
      VT(2) = T/4.0D0
      IF(VT(2).GT.XMAXM) VT(2)=XMAXM
      T = VT(2)
      DO 20 I=1,N
   20 X(I) = X(I)+T*P(I)*CONST(I)*SCALE
      CALL COMPFG(X,.TRUE.,ESCF,.TRUE.,EFS,.TRUE.)
      PHI(2)=DOT(EFS,EFS,N)
      CALL EXCHNG(PHI(2),SQSTOR,ENERGY,ESTOR,X,XSTOR,T,ALFS,N)
      DO 30 I=1,M
   30 GSTOR(I)=EFS(I)
      IF (PHI(1) .LE. PHI(2)) THEN
         VT(3) = -VT(2)
         LEFT = 3
         CENTER = 1
         RIGHT = 2
      ELSE
         VT(3)=2.0D0*VT(2)
         LEFT = 1
         CENTER = 2
         RIGHT = 3
      ENDIF
      TLAST = VT(3)
      T = TLAST-T
      DO 40 I=1,N
   40 X(I) = X(I)+T*P(I)*CONST(I)*SCALE
      FLAST=PHI(2)
      CALL COMPFG(X,.TRUE.,ESCF,.TRUE.,EFS,.TRUE.)
      F=DOT(EFS,EFS,N)
      IF(F.LT.SQSTOR) CALL EXCHNG(F,SQSTOR,ENERGY,ESTOR,X,
     1XSTOR,T,ALFS,N)
      DO 50 I=1,M
   50 GSTOR(I)=EFS(I)
      IF(F.LT.FIN) LOWER = .TRUE.
      NCOUNT = NCOUNT+2
      PHI(3) = F
      IF (IPRINT) 60,70,70
   60 WRITE (6,310) VT(1),PHI(1),VT(2),PHI(2),VT(3),PHI(3)
   70 MXCT=MXCNT2
      DO 250 ICTR=3,MXCT
         XMAXM=XMAXM*3.D0
         ALPHA = VT(2) - VT(3)
         BETA = VT(3) - VT(1)
         GAMMA = VT(1)-VT(2)
         IF(ALPHA.EQ.0.D0)ALPHA=1.D-20
         IF(BETA.EQ.0.D0)BETA=1.D-20
         IF(GAMMA.EQ.0.D0)GAMMA=1.D-20
         ABG =-(PHI(1)*ALPHA+PHI(2)*BETA+PHI(3)*GAMMA)/ALPHA
         ABG=ABG/BETA
         ABG=ABG/GAMMA
         ALPHA=ABG
         BETA = ((PHI(1)-PHI(2))/GAMMA)-ALPHA*(VT(1)+VT(2))
         IF (ALPHA)  80,80,110
   80    IF (PHI(RIGHT) .GT. PHI(LEFT))  GO TO 90
         T = 3.0D0*VT(RIGHT)-2.0D0*VT(CENTER)
         GO TO 100
   90    T = 3.0D0*VT(LEFT)-2.0D0*VT(CENTER)
  100    S=T-TLAST
         T=S+TLAST
         GO TO 150
  110    T = -BETA/(2.0D0*ALPHA)
         S=T-TLAST
         IF (S) 120,260,130
  120    AMDIS=VT(LEFT)-TLAST-XMAXM
         GO TO 140
  130    AMDIS=VT(RIGHT)-TLAST+XMAXM
  140    IF(ABS(S).GT.ABS(AMDIS)) S=AMDIS
         T=S+TLAST
  150    CONTINUE
         IF(ICTR.GT.3.AND.ABS(S*XMINM).LT.XCRIT) THEN
            IF( DEBUG )
     1    WRITE(6,'('' EXIT DUE TO SMALL PROJECTED STEP'')')
            GO TO 260
         ENDIF
         T=S+TLAST
         DO 160 I=1,N
  160    X(I) = X(I)+S*P(I)*CONST(I)*SCALE
         FLAST=F
         CALL COMPFG(X,.TRUE.,ESCF,.TRUE.,EFS,.TRUE.)
         F=DOT(EFS,EFS,N)
         IF(F.LT.SQSTOR) CALL EXCHNG(F,SQSTOR,ENERGY,ESTOR,X,XSTOR,
     1T,ALFS,N)
         DO 170 I=1,M
  170    GSTOR(I)=EFS(I)
         IF(F.LT.FIN) LOWER = .TRUE.
         NCOUNT = NCOUNT+1
         IF (IPRINT) 180,190,190
  180    WRITE (6,320) VT(LEFT),PHI(LEFT),VT(CENTER),PHI(CENTER),
     1  VT(RIGHT),PHI(RIGHT),T,F
  190    CONTINUE
C
C    TEST FOR EXCITED STATES AND POTHOLES
C
         ITRAP=0
         IF(ABS(VT(CENTER)).GT.1.D-10) GOTO 200
         IF(ABS(T)/(ABS(VT(LEFT))+1.D-15).GT.0.3333) GOTO 200
         IF(2.5D0*F-PHI(RIGHT)-PHI(LEFT).LT.0.5D0*PHI(CENTER)) GOTO 200
C
C   WE ARE STUCK ON A FALSE MINIMUM
C
         ITRAP=1
         GOTO 260
  200    CONTINUE
*
* NOW FOR THE MAIN STOPPING TESTS.  LOCMIN WILL STOP IF:-
*     THE ERROR FUNCTION HAS BEEN REDUCED, AND
*     THE RATE OF DROP OF THE ERROR FUNCTION IS LESS THAN 0.5% PER STEP
*     AND
*     (A) THE RATIO OF THE PROPOSED STEP TO THE TOTAL STEP IS LESS THAN
*         EPS,   OR
*     (B) THE LAST DROP IN ERROR FUNCTION WAS LESS THAN 5%OFTHETOTALDROP
*         DURING THIS CALL TO LOCMIN.
*
         IF(DEBUG)WRITE(6,'('' F/FLAST'',F13.6)')F/FLAST
         IF( LOWER  .AND. F/FLAST .GT. 0.995D0) THEN
            IF((ABS(T-TLAST).LE.EPS*ABS(T+TLAST)+TEE)) THEN
               IF( DEBUG )
     1       WRITE(6,'('' EXIT AS STEP IS ABSOLUTELY SMALL '')')
               GO TO 260
            ENDIF
            SUM=MIN(ABS(F-PHI(1)),ABS(F-PHI(2)),ABS(F-PHI(3)))
            SUM2=(FIN-SQSTOR)*0.05D0
            IF(SUM .LT. SUM2) THEN
               IF( DEBUG )
     1        WRITE(6,'('' EXIT DUE TO HAVING REACHED BOTTOM'')')
               GOTO 260
            ENDIF
         ENDIF
         TLAST = T
         IF ((T .GT. VT(RIGHT)) .OR. (T .GT. VT(CENTER) .AND. F .LT.
     1  PHI(CENTER)) .OR. (T .GT. VT(LEFT) .AND. T .LT. VT(CENTER) .AND.
     2  F .GT. PHI(CENTER)))  GO TO 210
         VT(RIGHT) = T
         PHI(RIGHT) = F
         GO TO 220
  210    VT(LEFT) = T
         PHI(LEFT) = F
  220    IF (VT(CENTER) .LT. VT(RIGHT))  GO TO 230
         I = CENTER
         CENTER = RIGHT
         RIGHT = I
  230    IF (VT(LEFT) .LT. VT(CENTER))  GO TO 240
         I = LEFT
         LEFT = CENTER
         CENTER = I
  240    IF (VT(CENTER) .LT. VT(RIGHT))  GO TO 250
         I = CENTER
         CENTER = RIGHT
         RIGHT = I
  250 CONTINUE
  260 CONTINUE
      CALL EXCHNG(SQSTOR,F,ESTOR,ENERGY,XSTOR,X,ALFS,T,N)
      DO 270 I=1,M
  270 EFS(I)=GSTOR(I)
      SSQ=(F)
      ALF=T
      IF (T) 280,300,300
  280 T = -T
      DO 290 I=1,N
  290 P(I) = -P(I)
  300 CONTINUE
      ALF=T
      RETURN
  310 FORMAT(' ---LOCMIN'/5X,'LEFT   ...',2F19.6/5X,'CENTER ...',
     1  2F19.6/5X,'RIGHT  ...',2F19.6/' ')
  320 FORMAT(5X,'LEFT   ...',2F19.6/5X,'CENTER ...',2F19.6/5X,
     1  'RIGHT  ...',2F19.6/5X,'NEW    ...',2F19.6/' ')
      END

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
Papa Johns was one of the better inexpensive pizza joints in Davis. Theyre a part of a national chain that lags only behind Little Caesers and Pizza Hut Express Pizza Hut in size. Not to be confused with Papa Murphys. Couponsaplenty arrive in the mail, making ordering here an even cheaper proposition. You can order online from their http://www.papajohns.com website (though you can now use some selected print coupons if you do, as well as new offers on the website.)  please note online ordering stops at 10:30pm now, so any order placed after that should not show up at the restaurant.

They had dipping sauces (garlic butter, pizza sauce) which are included with the pizzas as well as Bread Sticks, they also have (cheese sauce, buffalo sauce, blue cheese sauce, honey mustard sauce, and barbecue sauce), and they include a wiki:wikipedia:Pepperoncini with a pizza.

They also have wings and chicken strips which can come unflavored or with barbecue, buffalo, and honey chipotle sauces as well as 2 dipping sauce cups.

They have added small personal salads to their menu for $3.99, with the option of Italian or Ranch dressing.

They are now featuring a new improved Cinnapie which is the size of a small pizza 10 for $3.99.

They NOW charge $1.25 for delivery.  They just lowered the delivery charge just for Davis!

On Sunday, March 4th, 2007, they offered in a coupon given to past customers a large 1 topping pizza for $5.99, pickup only for Customer Appreciation Day, and as many pizzas as wanted could be had at that price.
  June 6, 2009 as well TusharRawat TR

Their 200910 junk mail ad. has a coupon for a $5.99 Medium 1Topping pizza but it says online only, so if you dont have FTP (Food Transfer Protocol), youre SOL (Outta Luck).
 Actually its not a carryout only, you place the order online but it is for delivery!! So you are IN luck!  Users/PapaJ
   The point is, not everyone has access to the internet at all hours, believe it or not. If its 9 PM, youre hungry, and cant get online, youre not going to be able to use the coupon.  Users/WilliamLewis

They have a current Carryout promotion for an Extra Large 1Topping pizza for $8.88.

20050116 03:14:46 nbsp I used to work at a Papa Johns  theyre good about food quality, but a little pricy unless you use a coupon. I suggest trying out the cheese sticks, they were my favorite. Users/GiladGurantz

20050116 12:46:10 nbsp the best coupon is the one in the green coupon machine book that you get in the mail Users/RohiniJasavala

20050120 15:32:32 nbsp Currently running a special for their 21st Anniversay: Buy a Large Combo Type Pizza ( I suggest the Hawaiian BBQ! ) get a Large 1 Topping for Free.  2 Large Pizzas for under $20 Delivered ($18.85 or something like that) Users/TarZxf

20050219 22:21:52 nbsp Theyll deliver to campus, my friends and I have had them deliver to rooms in Wellman Hall Wellman and Olson Hall Olson. Users/RyanCastellucci

20050313 14:01:04 nbsp Ive gotten this a few times and have been consistently dissapointed by the dough. Users/ArlenAbraham

20050702 21:13:33 nbsp In Neal Stephensons epic wiki:WikiPedia:Cyberpunk cyberpunk novel wiki:WikiPedia:Snow_Crash Snow Crash, the United States  now in shambles  can only do two things better than any other nation:  Deliever pizza and write microcode.  Those who work in the pizza industry attend a specific University for four years, and if a pizza isnt delievered in 30 minutes or less, news is flashed to the Pizza Headquarters.  Immediately, the owner of the corporation calls the customer within five minutes, apologizes profusely, flies to the customers home in a jet helicopter, apologizes more, and gives him a free trip to Italy.

If Papa Johns is any indication, we are unreasonably far from this fastpaced cyberpunk reality.  My pizza, which I ordered online at 11PM, was never delivered.  It turns out that you can order pizza through their website up until 11:30PM, but their store closes at 11PM.  Thus, all intermediate orders get lost.  I called them the next day, starving, and told them I never got my pizza and that my credit card had been charged.  They didnt apologize and asked me what kind of pizza I had ordered.  Cheese, I remarked.  In about 30 minutes they came with my pizza  13 hours late  along with a copy of the bill.

Also, their garlic sauce is way too buttery. Users/PhilipNeustrom

20050703 04:17:08 nbsp I hate online orders. They just show up, so if youre not looking at the screen constantly, they get totally ignored. Just call!  Note from a soontobemanager at Papa Johns in San Jose. Users/GiladGurantz

20050703 13:25:25 nbsp Yeah, I find online orders are a terrible way to go, too.  NOBODY does online orders well yet.  In fact, if I have to make an online order to get a special discount, Ill usually follow up with a phone call to make sure they got it. Users/JaimeRaba
   If Papa Johns offers the service they should honor it.  If online orders are stupid and the website allows you to order pizza past the closing time of the store, then maybe they should stop offering it!  I only did the online order because I wasnt sure you could pay with credit card when you get it delievered.  Can you?
of course you can, if they dont have a portable swiper theyll take your info over the phone and punch it in by hand.

20050825 21:55:24 nbsp I order pizza online and it was easy and painless. The pizza showed up on time and so dont badmouth the online service. Give it a try. Users/SimonFung

20050918 00:25:31 nbsp Since we now live 5 minutes walking distance from here, I decided to try Papa Johns. The pizza was OK, but Users/BlancheNonken and I felt the sauce was too sweet. I consider that a diet issue if theyre sweetening it in some way (adding refined sugar or some such). Users/JeffreyNonken

20051001 15:16:44 nbsp The packeted parmesean cheese is slightly more bitter compared to other pizza chains. Users/AlexanderHo

20051112 00:56:12 nbsp Best pizza joint in town! If ever there was a question about a pizza to get, Papa Johns here we come. Doughy wonderful crust, good amount of sauce, perfect balance of cheese and other toppings. Much higher quality pizza than that Roundtable and Dominos stuff they try to pass off as pizza.  Users/TarynHornen

20060214 15:10:18 nbsp The website ordering isnt well thought out by the corporation, and is a giant pain in the ass for those who work there. Calling to make sure the order was recieved after you place it is a really good idea. Users/DaveCar

20061104 14:37:51 nbsp Round Table has higher quality pizza and service than Papa Johns. Users/TusharRawat

20061108 21:00:10 nbsp they have a buy one get one free discount on tuesday. i like their bbq chicken pizza and dont forget the garlic sauce Users/AnnThiNguyen

20061210 18:18:43 nbsp Everyone is entitled to their own opinion Tushar. If your claims are fact, why dont you provide some backup data?  Users/DaveCar
  That was an opinion. Notice that I did not say the pizza was bad. I indicated that I like RTP better.

20070116 23:38:55 nbsp I just had a pizza from here and it tasted like something Id buy at safeway.  That is to say it was way too doughy, and the cheese was rubbery.  Oh how I long for the Pizza Azzuro of home.  Why isnt there a decent traditional Italian style pizza place around here? Users/NicholasKnoblauch

20070127 22:57:57 nbsp i saw something about them on the food network and will no longer purchase pizza any where else. truly inspiring :O) on a side note i also heard that the teenage mutant ninja turtles purchase their pizza from these guys too :O) Users/EmilyTung

20070425 16:50:43 nbsp We all like Papa Johns pizza, the kids working there are always helpful and courteous, even when the customers are being NASTY!  Especially one kid working there, I think his name is Adam, he always keeps his cool and is super polite too  a nice touch in a fast food place like this. Users/DavisGradStudent

20071103 17:20:59 nbsp They have a $7.99 for Large and $9.99 for extra large onetopping if you go pick it up. I think its a daily offer for the time being. Users/TusharRawat

20071105 13:18:31 nbsp I just love the pizza so much and the price is not bad at all. But, I was very dissapointed with the employees, they were not really friendly. I said thanks to them and say good night, and they did not say anything, they just ignored me. I believe, I did nothing wrong, I felt so bad because I just told by myself. Can you guys smile a little bit please?  Users/viviancute

20080123 21:22:30 nbsp Vomitrocious. Users/thelonepiper555

20080213 20:22:05 nbsp I have been working there since October of 2007. We used to be ok but now we are under a terrible new manager. none of the managers are skilled enough to do their jobs properly (Except Maria) and one of them is a god damned illiterate. Illiterate and still makes $1.25 more than me. He only got the job because he knows the manager through the managers cousin. For the love of humanity, spend your hard earned money elsewhere. Users/AndresRuiz

20080719 22:22:19 nbsp I just got home from a 23 mile bike ride and was hungry for some pizza. Looked through the assortment of pizza coupons I had and decided to try Papa Johns. Called in my order, was put on hold, then some guy (i forgot his name) took my order. Couldnt get the coupon right, but eventually my order was placed and I was told it would be ready in 15 min. Went to go pick up the pizza and when I got there it wasnt made and they said it wasnt in the system They asked me which number I called, I couldnt remember the number off hand because it was on the coupon sheet and I called it from my house phone. They offered to make it for me but didnt offer any further discount so I decided to go elsewhere. Sure enough when I got home and checked the number I dialed it was right. I probably should have gone to Pizza Guys, Ive never had a problem with my order there.    Users/DavidHolcomb

20090117 21:40:10 nbsp This is the best pizza of the ridiculously cheap pizzas, as long as you have a coupon. They are also very speedy on the delivery (and as I remember, the delivery guys arent bad looking either). Its a regular greasy tradition in my house to buy cheesy bread  always a good decision. Users/caitlinpolly

20091030 15:34:33 nbsp There has been several changes over the years. With the current Executive Management our goal is to bring back the good days of pizza and strive for excellent customer service going forward. We appreciate all of the folks in Davis, especially the Students and the UCD Campus. Also we very much appreciate the folks who live in the local community year round. We hope to bring you a quality pizza with a good value. Users/PapaJ

20091103 20:36:43 nbsp Tha same pizza that was $7 last week is now $17. The cheapness used to be enough to overlook the crap sugary sauce. The only decent pizza in this town is Symposium and Village Pizza. Users/JerseyCity

20091130 07:43:01 nbsp they are good, they work hard and they are much cleaner than some other pizza places in town. people stop complaining about prices, yes a large pepperoni pizza is $15, but go ahead and get one from some other places with same quality then camper the price. and Andres i remember you working there too, you werent a good worker like some of our shift leaders, however you sucked! and also Maria sucks too as well, but sadly she can not go anywhere else. Users/bigO

20091205 16:11:58 nbsp For some reason in the last month the website for the Davis location does not accept coupon codes that I get in my email each week.  I called in and management apparently hasnt changed...they honored it over the phone for me instead. Users/OliviaY

20100325 15:44:01 nbsp Papa Johns delivery charge is now Only $1.25, and these great specials: Any Medium Pizza for Only $8.00, Any Large Pizza for Only $10, or Any Extra Large Pizza for Only $12. With prices like that, why not get two pizzas? Users/PapaJ

20100504 20:34:59 nbsp Ill tell you a reason not to get 2 pizzas... i called in just now to order a pizza for pick up, and they wouldnt take my order claiming they were too busy!  Ive had small problems multiple times from this papa johns, and theyre receiving less and less of my business.  Ive never heard of a pizza place doing this!  I like their pizza.... but come on... what place is too busy to take an order, especially when im asking to pick it up and not putting a drain on their delivery drivers.   Users/marshvz

20100525 22:57:22 nbsp Theres a sign posted on their door that says this location is now closed.  As in, closed for good...  Users/JaneGreene

20100526 18:01:22 nbsp It was dark when I rode by last night, but the signage was still up

man that shopping center is struggles on the side tenants, at least it is anchored with the coop Users/StevenDaubert

20100723 14:28:29 nbsp Say it aint so!! How else am I gonna get my Garlic Butter fix? Users/TrishAng

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
subroutine binmodels_py(wv1,wv2,dw,snpt,starmodel_wv,starmodel_flux,ld_coeff, &
	pnpt,planetmodel_wv,planetmodel_rprs, &
	bmax,bin_starmodel_wv,bin_starmodel_flux,bin_ld_coeff,bin_planetmodel_wv, &
	bin_planetmodel_rprs)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
use precision
implicit none
integer :: snpt,pnpt,bmax
real(double) :: wv1,wv2,dw

real(double) :: starmodel_wv(snpt),starmodel_flux(snpt),ld_coeff(snpt,4)
real(double) :: bin_starmodel_wv(bmax),bin_starmodel_flux(bmax),bin_ld_coeff(bmax,4)

real(double) :: planetmodel_wv(pnpt),planetmodel_rprs(pnpt)
real(double) :: bin_planetmodel_wv(bmax),bin_planetmodel_rprs(bmax)

integer :: i,b
real(double) :: s_wv
integer, dimension(:), allocatable :: bin_count

!CCCCCCCCCCCCCCCCCCCCCC
! Start with Star model
!CCCCCCCCCCCCCCCCCCCCCC

!Initialize array
bin_starmodel_wv=0.0d0
bin_starmodel_flux=0.0d0
bin_ld_coeff=0.0d0

!allocate array to count samples in each bin
allocate(bin_count(bmax))
bin_count=0

do i=1,snpt

	s_wv=starmodel_wv(i) !current wavelength
	b=int((s_wv-wv1)/dw)+1 !bin number

	if ((b.gt.0).and.(b.le.bmax))then
		bin_starmodel_wv(b)=bin_starmodel_wv(b)+starmodel_wv(i)
		bin_starmodel_flux(b)=bin_starmodel_flux(b)+starmodel_flux(i)
		bin_ld_coeff(b,:)=bin_ld_coeff(b,:)+ld_coeff(i,:)
		bin_count(b)=bin_count(b)+1
	endif

enddo

do b=1,bmax
	if (bin_count(b).gt.0) then
		bin_starmodel_wv(b)=bin_starmodel_wv(b)/dble(bin_count(b))
		bin_starmodel_flux(b)=bin_starmodel_flux(b)/dble(bin_count(b))
		bin_ld_coeff(b,:)=bin_ld_coeff(b,:)/dble(bin_count(b))
	else
		bin_starmodel_wv(b)=dw*(dble(b)+0.5)+wv1
	endif
enddo

!CCCCCCCCCCCCCCCCCCCCCC
! Bin Planet Model
!CCCCCCCCCCCCCCCCCCCCCC

bin_planetmodel_wv=0.0d0
bin_planetmodel_rprs=0.0d0

!reset bin_count
bin_count=0

do i=1,pnpt

    s_wv=planetmodel_wv(i) !current wavelength
	b=int((s_wv-wv1)/dw)+1 !bin number

	if ((b.gt.0).and.(b.le.bmax))then
		bin_planetmodel_wv(b)=bin_planetmodel_wv(b)+planetmodel_wv(i)
		bin_planetmodel_rprs(b)=bin_planetmodel_rprs(b)+planetmodel_rprs(i)
		bin_count(b)=bin_count(b)+1
	endif

enddo

do b=1,bmax
	if (bin_count(b).gt.0) then
		bin_planetmodel_wv(b)=bin_planetmodel_wv(b)/dble(bin_count(b))
		bin_planetmodel_rprs(b)=bin_planetmodel_rprs(b)/dble(bin_count(b))
	else
		bin_planetmodel_wv(b)=dw*(dble(b)+0.5)+wv1
	endif
enddo

return
end subroutine binmodels_py

max_stars_repo_name: AntoineDarveau/jwst-mtl

max_issues_repo_name: AntoineDarveau/jwst-mtl

max_forks_repo_name: AntoineDarveau/jwst-mtl


text:
Apartment Name is B Street Apartments.

Complex has fifteen 2bedroom townhouse style apartments, one 3bedroom apartment and one 1bedroom apartment. The complex only has hours set up by the Property Manager who lives in the complex at Apartment #1 and whose hours are from 36pm Mondays  Fridays and any other time by appointment.

 Management Notices

##YYYYMMDD  Mgmt can leave notes timely notes here about specials, discounts, promotions, announcements etc etc.  Do NOT post huge, glossy flyers.  This is not an advertisement space, it is a place to share timely information to prospective and current residents.

 Management can leave timely notes here.  Please dont spam!

 Features & Amenities

 Mobility
Public Transportation  Apartment is located on Unitrans E bus line.  Several Yolobus routes stop nearby at Russell or at the Memorial Union.  That Street and will take you to Sacramento, West Sacramento, Woodland and Sacramento International Airport

Bicycle  You can get to The Davis Bike Loop by going through campus.  The UC Davis campus is less than 2 block away.  There are bike lanes/connectors on all/some/these of the surrounding streets.

Pedestrian  For recreational adventures, walk to downtown.  For grocery, shopping, some services and dining, the CoOp Grocery Store is on G street. Farmers market is across the street in the park.

Parking  No guests are allowed to park within the complex.  That is strictly for residents.  You will be towed if you over extend your stay in the Black Bear Dining Parking Lot.  Street parking is available but limited.  Also bike parking is limited.

ADA Accessibility  The bottom floor units are ADA Accessible.

 Floor Plans
If you know the different apartment sizes with their corresponding rent, you can put the information in the following table for easy reference:

  To learn more about rental housing in Davis, check out our Housing Guide as well as our Apartments pages.

## uncomment the next line when it is time to create a comment archive
##Older comments can be found in the /Comment Archive comment archive.

Im sure people have perfectly fine apartments in this complex.  However, my apartment has termites and mold, which is apparently common in all of the apartments.  My garbage disposal broke, my microwave doesnt rotate and my light switches (the electrical wiring of the house) are on the fritz.  I talked to them two and a half weeks ago and they have yet to do anything about the termites, except bring in a maintenance man to look at the pile of dead carcasses on my window sills, floors and dishes.  Also he helped himself to my toilet and forgot to put the seat down.  Professional, much?  The property doesnt allow pets but the property manager has two dogs living with her.  The Management (Brittain Property Management) has received terrible reviews, which can be viewed on several different sites, including the one that is posted under Website at the top of the page, and I think Ill let those reviews speak for themselves.  Apparently they do not have an email, which I find sketchy as far as businesses go.  How does a business not have an email in the 21st Century? OR a website?? Parking is bad for visitors because Black Bear will tow if you stay long enough and youd need to buy an extra permit (one time fee of $12) from the city in order to park on the street.  Overall, this place is next to squalor.  The structural integrity, the management, the maintenance, and the services are terrible and the rent is ridiculously high.  The only thing positive I could say about this is the location, but Im sure you could move on over to A Street Apartments or Villa Verde instead.  If a Building Inspector came into this place (let alone a Health and Safety Inspector), this property would be foreclosed or, at the very least, the business would lose most of their credentials.  THAT is how bad it is.  Keep looking.

20120111 10:53:25 nbsp Really? I live here and the place is not as bad as the previous reviewer stated. Parking is tough for visitors and they really should have a webpage or email. But my place is pretty nice and the management is really quick in sending in repair requests. He is right about location. 1min bike ride to campus. Users/shalimar


text:
      double complex Hbit,Bbit
      common/WWbits/Hbit,Bbit
!$omp threadprivate(/WWbits/)

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
! { dg-do run }
! Tests the fix for pr28167, in which character array constructors
! with an implied do loop would cause an ICE, when used as actual
! arguments.
!
! Based on the testscase by Harald Anlauf  <anlauf@gmx.de>
!
  character(4), dimension(4) :: c1, c2
  integer m
  m = 4
! Test the original problem
  call foo ((/( 'abcd',i=1,m )/), c2)
  if (any(c2(:) .ne. (/'abcd','abcd', &
                            'abcd','abcd'/))) call abort ()

! Now get a bit smarter
  call foo ((/"abcd", "efgh", "ijkl", "mnop"/), c1) ! worked previously
  call foo ((/(c1(i), i = m,1,-1)/), c2)            ! was broken
  if (any(c2(4:1:-1) .ne. c1)) call abort ()

! gfc_todo: Not Implemented: complex character array constructors
  call foo ((/(c1(i)(i/2+1:i/2+2), i = 1,4)/), c2)  ! Ha! take that..!
  if (any (c2 .ne. (/"ab  ","fg  ","jk  ","op  "/))) call abort ()

! Check functions in the constructor
  call foo ((/(achar(64+i)//achar(68+i)//achar(72+i)// &
               achar(76+i),i=1,4 )/), c1)           ! was broken
  if (any (c1 .ne. (/"AEIM","BFJN","CGKO","DHLP"/))) call abort ()
contains
  subroutine foo (chr1, chr2)
    character(*), dimension(:) :: chr1, chr2
    chr2 = chr1
  end subroutine foo
end

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


text:
      SUBROUTINE newt(x,n,check)
      INTEGER n,nn,NP,MAXITS
      LOGICAL check
      REAL x(n),fvec,TOLF,TOLMIN,TOLX,STPMX
      PARAMETER (NP=40,MAXITS=200,TOLF=1.e-4,TOLMIN=1.e-6,TOLX=1.e-7,
     *STPMX=100.)
      COMMON /newtv/ fvec(NP),nn
      SAVE /newtv/
CU    USES fdjac,fmin,lnsrch,lubksb,ludcmp
      INTEGER i,its,j,indx(NP)
      REAL d,den,f,fold,stpmax,sum,temp,test,fjac(NP,NP),g(NP),p(NP),
     *xold(NP),fmin
      EXTERNAL fmin
      nn=n
      f=fmin(x)
      test=0.
      do 11 i=1,n
        if(abs(fvec(i)).gt.test)test=abs(fvec(i))
11    continue
      if(test.lt..01*TOLF)then
        check=.false.
        return
      endif
      sum=0.
      do 12 i=1,n
        sum=sum+x(i)**2
12    continue
      stpmax=STPMX*max(sqrt(sum),float(n))
      do 21 its=1,MAXITS
        call fdjac(n,x,fvec,NP,fjac)
        do 14 i=1,n
          sum=0.
          do 13 j=1,n
            sum=sum+fjac(j,i)*fvec(j)
13        continue
          g(i)=sum
14      continue
        do 15 i=1,n
          xold(i)=x(i)
15      continue
        fold=f
        do 16 i=1,n
          p(i)=-fvec(i)
16      continue
        call ludcmp(fjac,n,NP,indx,d)
        call lubksb(fjac,n,NP,indx,p)
        call lnsrch(n,xold,fold,g,p,x,f,stpmax,check,fmin)
        test=0.
        do 17 i=1,n
          if(abs(fvec(i)).gt.test)test=abs(fvec(i))
17      continue
        if(test.lt.TOLF)then
          check=.false.
          return
        endif
        if(check)then
          test=0.
          den=max(f,.5*n)
          do 18 i=1,n
            temp=abs(g(i))*max(abs(x(i)),1.)/den
            if(temp.gt.test)test=temp
18        continue
          if(test.lt.TOLMIN)then
            check=.true.
          else
            check=.false.
          endif
          return
        endif
        test=0.
        do 19 i=1,n
          temp=(abs(x(i)-xold(i)))/max(abs(x(i)),1.)
          if(temp.gt.test)test=temp
19      continue
        if(test.lt.TOLX)return
21    continue
      pause 'MAXITS exceeded in newt'
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/newt.for


text:
\ Chad Tiny Encryption Algorithm							9/19/20 BNE

\ XTEA-inspired crypto algorithm for Chad.
\ Uses 2-cell data and a 4-cell key.
\ For 18-bit cells, it's 36-bit data and a 72-bit key.

\ CTEA uses 71 instructions of code space.
\ The numeric conversion pointer `hld` is used as temporary storage.

there

\ Note: Key lengths over 56-bit may be subject to export controls.
\ The last cell in the table is 0 to make it a 54-bit key (if 18-bit cell).

align 4 cells buffer: CTkey
$179B9 equ CTdelta
18 equ CTrounds \ about 110 cycles per round

: CTshift  ( v -- v v' sum sum )
    dup 2* 2* 2* 2*  over 2/ 2/ 2/ 2/ 2/  xor  \ could be custom instruction
	over +	hld @  dup
;
: CTcalcA  CTshift [ ;
: CTdokey  3 and cells CTkey + @  +  xor  rot swap ;
: CTcalcB  CTshift  swapb 2/ 2/ 2/  CTdokey ;

: encipher  ( v0 v1 -- v0' v1' )
    0 hld !
    CTrounds for
		CTcalcA +
		CTdelta hld +!
		CTcalcB +
	next
	swap
;
: decipher  ( v0 v1 -- v0' v1' )
    [ CTdelta CTrounds * ] literal hld !
    CTrounds for
	    CTcalcB -
		[ CTdelta negate ] literal  hld +!
	    CTcalcA -
	next
;

there swap - . .( instructions used by ctea) cr

max_stars_repo_name: uho/chad

max_issues_repo_name: uho/chad

max_forks_repo_name: uho/chad


text:
program enkf_main
!$$$  main program documentation block
!
! program:  enkf_main                  high level driver program for
!                                      ensemble kalman filter (EnKF).
!
! prgmmr: whitaker         org: esrl/psd               date: 2009-02-23
!
! abstract: This is the main program for the EnKF code. It does the following:
!           a) initialize MPI, read namelist from enkf.nml on each task.
!           b) reads observations, observation priors and associated
!              metadata for each ensemble member (from diag file
!              output by GSI forward operator code). Print innovation
!              statistics for prior.
!           c) read horizontal grid information (lat/lon of each grid point) and
!              pressure at each grid point/vertical level.
!           d) decomposition of horizontal grid points and observation
!              priors to minimize load imbalance.
!           e) read state variables for each ensemble member, distribute
!              to each task.
!           f) update state variables, observation priors and radiance
!              bias correction coefficients with EnKF. The observation
!              prior perturbations updated by each task are read in
!              from a temp file created in step b).
!           g) inflate posterior ensemble perturbations proportional
!              to the amount that the variance was reduced by the analysis.
!           h) write out the analysis ensemble and updated radiance bias
!              correction coefficients. Print innovation statistics for
!              posterior.
!           i) deallocate all allocatable arrays, finalize MPI.
!
! program history log:
!   2009-02-23  Initial version.
!   2011-06-03  Added the option for LETKF.
!   2016-02-01  Initialize mpi communicator for IO tasks (1st nanals tasks).
!   2017-05-20  Allows for calc of prior and posterior residual stats to
!               include only obs actually assimilated.
!
! usage:
!   input files:
!     sfg_YYYYMMDDHH_fhr06_mem* - first guess ensemble members, plus
!                     ensemble mean (sfg_YYYYMMDDHH_fhr06_ensmean).
!     satbias_angle - satellite angle dependent file
!     satbias_in    - satellite bias correction coefficient file
!     satinfo       - satellite channel info file
!     convinfo      - convential data (prepufr) info file
!     ozinfo        - ozone retrieval info file
!     diag_YYYYMMDDHH_ges_mem*  - observation diagnostic files for each ensemble member
!                     created GSI forward operator.
!     hybens_locinfo - if parameter readin_localization is true, contains
!                      vertical profile of horizontal and vertical localization
!                      length scales.
!
!   output files:
!     sanl_YYYYMMDDHH_mem* - analysis ensemble members. A separate program
!                            may be run to add system noise to these files.
!     covinflate.dat - multiplicative inflation (inflate_ens in module inflation).
!     satbias_out    - output satellite bias correction file.
!
! comments:
!
! This program is run after the forward operator code is run on each ensemble
! member to create the diag*mem* input files.
!
! attributes:
!   language: f95
!
!$$$

 use kinds, only: r_kind,r_double,i_kind
 ! reads namelist parameters.
 use params, only : read_namelist,letkf_flag,readin_localization,lupd_satbiasc,&
                    numiter, nanals, lupd_obspace_serial,&
                    fso_flag, fso_cycling, stats_usedob_only
 ! mpi functions and variables.
 use mpisetup, only:  mpi_initialize, mpi_initialize_io, mpi_cleanup, nproc, &
                      numproc, mpi_wtime
 ! obs and ob priors, associated metadata.
 use enkf_obsmod, only : readobs, obfit_prior, obsprd_prior, &
                    deltapredx, nobs_sat, obfit_post, obsprd_post, &
                    obsmod_cleanup, biasprednorminv
 ! innovation statistics.
 use innovstats, only: print_innovstats
 ! grid information
 use gridinfo, only: getgridinfo, gridinfo_cleanup, npts,lonsgrd,latsgrd
 ! model state vector
 use statevec, only: read_ensemble, write_ensemble, statevec_cleanup
 ! load balancing
 use loadbal, only: load_balance, loadbal_cleanup
 ! enkf update
 use enkf, only: enkf_update
 ! letkf update
 use letkf, only: letkf_update
 ! radiance bias correction coefficients.
 use radinfo, only: radinfo_write, predx, jpch_rad, npred
 ! FSO update
 use fso, only: fso_update
 ! Observation sensitivity usage
 use enkf_obs_sensitivity, only: init_ob_sens, print_ob_sens, destroy_ob_sens, &
                                 loc_advection, read_ob_sens
 ! posterior ensemble inflation.
 use inflation, only: inflate_ens
 ! initialize radinfo variables
 use radinfo, only: init_rad, init_rad_vars
 use omp_lib, only: omp_get_max_threads

 implicit none
 integer(i_kind) j,n,nth
 real(r_double) t1,t2
 logical no_inflate_flag
 ! initialize MPI.
 call mpi_initialize()
 if (nproc==0) call w3tagb('ENKF_ANL',2011,0319,0055,'NP25')

 ! Initial radinfo variables (some flags may be over-ridden in enkf namelist)
 call init_rad()

 ! read namelist.
 call read_namelist()

 ! initialize MPI communicator for IO tasks.
 call mpi_initialize_io(nanals)

 ! Initialize derived radinfo variables
 call init_rad_vars()

 nth= omp_get_max_threads()
 if(nproc== 0)write(6,*) 'enkf_main:  number of threads ',nth

 ! read horizontal grid information and pressure fields from
 ! 6-h forecast ensemble mean file.
 call getgridinfo()

 if (.not. fso_flag) then
   ! read obs, initial screening.
   t1 = mpi_wtime()
   call readobs()
   t2 = mpi_wtime()
   if (nproc == 0) print *,'time in read_obs =',t2-t1,'on proc',nproc

   ! print innovation statistics for prior on root task.
   if (nproc == 0) then
     if(.not.stats_usedob_only) then
        print *,'innovation statistics for prior:'
        call print_innovstats(obfit_prior, obsprd_prior)
     endif
   end if
 else
   ! Read data for observation sensitivity computation
   call read_ob_sens()
 end if

 ! Initialization for writing
 ! observation sensitivity files
 if(fso_cycling) call init_ob_sens()

 ! read in vertical profile of horizontal and vertical localization length
 ! scales, set values for each ob.
 if (readin_localization) call read_locinfo()

 ! do load balancing (partitioning of grid points, observations among
 ! processors)
 t1 = mpi_wtime()
 call load_balance()
 t2 = mpi_wtime()
 if (nproc == 0) print *,'time in load_balance =',t2-t1,'on proc',nproc

 ! read in ensemble members, distribute pieces to each task.
 t1 = mpi_wtime()
 call read_ensemble()
 t2 = mpi_wtime()
 if (nproc == 0) print *,'time in read_ensemble =',t2-t1,'on proc',nproc

 ! advection of localization function
 if(fso_flag) then
    call init_ob_sens()
    t1 = mpi_wtime()
    call loc_advection()
    t2 = mpi_wtime()
    if (nproc == 0) print *,'time in loc_advection =',t2-t1,'on proc',nproc
 end if

 t1 = mpi_wtime()
 ! state and bias correction coefficient update iteration.
 if(fso_flag) then
    call fso_update()
 else if(letkf_flag) then
    ! do ob space update using serial filter if desired
    if (lupd_obspace_serial) call enkf_update()
    call letkf_update()
 else
    call enkf_update()
 end if
 t2 = mpi_wtime()
 if (nproc == 0) print *,'time in enkf_update =',t2-t1,'on proc',nproc

 ! Output non-inflated
 ! analyses for FSO
 if(fso_cycling) then
   no_inflate_flag=.true.
   t1 = mpi_wtime()
   call write_ensemble(no_inflate_flag)
   t2 = mpi_wtime()
   if (nproc == 0) print *,'time in write_ensemble wo/inflation =',t2-t1,'on proc',nproc
 end if

 ! posterior inflation.
 if (.not. fso_flag) then
   t1 = mpi_wtime()
   call inflate_ens()
   t2 = mpi_wtime()
   if (nproc == 0) print *,'time in inflate_ens =',t2-t1,'on proc',nproc
 end if

 ! print observation sensitivity summaries on root task.
 if(fso_cycling .or. fso_flag) call print_ob_sens()

 ! print innovation statistics for posterior on root task.
 if (.not. fso_flag .and. nproc == 0 .and. numiter > 0) then
    if (stats_usedob_only) then
       print *,'innovation statistics for prior:'
       call print_innovstats(obfit_prior, obsprd_prior)
    endif
    print *,'innovation statistics for posterior:'
    call print_innovstats(obfit_post, obsprd_post)
 ! write out bias coeffs on root.
    if (nobs_sat > 0 .and. lupd_satbiasc) then
       ! re-scale bias coefficients.
       do j=1,jpch_rad
           do n=1,npred
              predx(n,j) = predx(n,j)*biasprednorminv(n)
           enddo
       enddo
       call radinfo_write(nproc)
    end if
 end if

 ! free memory (radinfo memory freed in radinfo_write)
 ! and write out analysis ensemble.
 call obsmod_cleanup()

 if (.not. fso_flag) then
   no_inflate_flag=.false.
   t1 = mpi_wtime()
   call write_ensemble(no_inflate_flag)
   t2 = mpi_wtime()
   if (nproc == 0) print *,'time in write_ensemble w/inflation =',t2-t1,'on proc',nproc
 end if

 call gridinfo_cleanup()
 call statevec_cleanup()
 call loadbal_cleanup()
 call destroy_ob_sens()

 ! write log file (which script can check to verify completion).
 if (nproc .eq. 0) then
    call write_logfile()
 endif

 ! finalize MPI.
 if (nproc==0) call w3tage('ENKF_ANL')
 call mpi_cleanup()

end program enkf_main

max_stars_repo_name: GEOS-ESM/GEOSadas

max_issues_repo_name: GEOS-ESM/GEOSadas

max_forks_repo_name: GEOS-ESM/GEOSadas


text:
!***********************************************************************
!                                                                      *
      SUBROUTINE CONNECT
!                                                                      *
!   The position of an orbital in the merged list is connected to      *
!   the positions in the initial and final state lists                 *
!                                                                      *
!   Written by Per Jonsson                                             *
!                                                                      *
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  13:35:54   1/ 3/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE biorb_C
      USE orb_C
      IMPLICIT NONE
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: I, J
!-----------------------------------------------
!
!ww      INTEGER PNTRIQ
!     POINTER (PNTRIQ,RIQDUMMY)
!
!
!     Initialize
!
      NNII(:NW) = 0.D0
      NNFF(:NW) = 0.D0
!
!   Loop over the orbitals in the merged list
!
      DO I = 1, NW
         DO J = 1, NWII
            IF (NP(I)/=NPII(J) .OR. NAK(I)/=NAKII(J)) CYCLE
            NNII(I) = J
         END DO
         DO J = 1, NWFF
            IF (NP(I)/=NPFF(J) .OR. NAK(I)/=NAKFF(J)) CYCLE
            NNFF(I) = J
         END DO
      END DO

      RETURN
      END SUBROUTINE CONNECT

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
cc
cc    Compile: gfortran vvv_flx2mag_new.f -o vvv_flx2mag -Lcfitsio -lcfitsio -lm
cc    Compile on VVV server as: gfortran vvv_flx2mag.f -o vvv_flx2mag -L/sw/lib -lcfitsio -lm
cc
cc    This is a modified version of the code fitsio_cat_list written by CASU.
cc    This code computes the magnitudes and their errors in the first 7 apertures.
cc
cc    Written by Istvan Dekany
cc    idekany (a t) astro (d o t) puc (d o t) cl
cc
      parameter (ncmax=80,nrmax=2000000,napers=7)
      real*8 tpa,tpd,tand,secd,tpi,pi
      real*8 xx,yy,xi,xn
      real*8 flux(napers),fluxerr(napers),fixed(napers)
      real*4 rec(nrmax,ncmax)
      real*4 magzpt,magzrr,maglim
      real*8 magobj(napers),magerr(napers)
      real*4 bin(35,35)
      real*4 apcor(napers)
      real*4 delapcor(napers)
      character*80 infile,outfile
      character*80 argv,errmsg,dateobs,utstart,tablename
      character*72 comment
      character*68 value
      character*16 ttype(ncmax),tunit(ncmax),tform(ncmax)
      character*11 filter
      character*1 csign,czero
      integer status,blocksize
      logical first,illum
      common /illum/ first,bin
      common /rd/ tpa,tpd,tand,secd,tpi,a,b,c,d,e,f,projp1,projp3,
     1            projp5,iwcstype

      li=3
      lr=5
      lt=6
      pi=4.0*atan(1.0)
      tpi=8.0*atan(1.0)
      czero='0'
      first=.true.
      illum=.false.                         ! default is no illumination table

c *** try for presence of command line arguments
      m=iargc()
c *** if correct number of arguments try and interpret them
      if(m.eq.2.or.m.eq.3)then
        call getarg(1,argv)
        infile=argv
        call getarg(2,argv)
        outfile=argv
        if(m.eq.3)then
          call getarg(3,argv)
          tablename=argv
          illum=.true.
        endif
      else
        print *,' '
        print *,'Arguments: input output [illum. corr. table]'
        print *,' '
        stop
      endif
c *** open fits file and read in catalogue info
      iunit=1
      status=0
      call ftopen(iunit,infile,0,blocksize,status)        ! 0 readonly option
      if(status.ne.0)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
        ier=1
        status=0
        call ftclos(iunit,status)
      endif
      call ftthdu(iunit,nhdus,status)
      nmefs=max(1,nhdus-1)

c *** clobber output file if exists
      call clobber(3,outfile)
      open(unit=3,file=outfile,status='new')

c *** see if VIRCAM or VST and get some info from PHU
      iesocam=0
      call ftgkey(iunit,'INSTRUME',value,comment,status)
      if(status.eq.202)then
        status=0
      else
        if(index(value,'VIRCAM').gt.0)iesocam=1
        if(index(value,'OMEGACAM').gt.0)iesocam=2
        if(iesocam.ne.0)then
          call ftgkys(iunit,'ESO INS FILT1 NAME',filter,comment,
     1                status)
          call ftgkye(iunit,'ESO TEL AIRM END',airmass,comment,
     1                status)
          call ftgcrd(iunit,'DATE-OBS',dateobs,status)
          utstart='UTST    ='
c *** try and read the exptime from here as a precaution
          call ftgkye(iunit,'EXPTIME',expsave,comment,status)
          if(status.eq.202)then
            expsave=0.0
            status=0
          endif
        endif
      endif

c *** now go through each extension
      do mef=1,nmefs
      call fitsio_rtable(rec,iunit,mef,nrows,ncols,apcor,delapcor,
     1  apcorpk,flim,filter,exptime,utstart,dateobs,percorr,magzpt,
     2  magzrr,airmass,extinct,skyvar,gain,saturate,pltscl,theta,
     3  iesocam)

c *** check for sensible exptime value
      if(exptime.eq.0.0)then
        exptime=expsave                                   ! see if PHU value ok
        if(exptime.eq.0.0)then
          print *,'exposure time information missing'
          stop
        endif
      endif

c *** some magnitude info
      magzpt=magzpt-(airmass-1.0)*extinct+2.5*log10(exptime)
      maglim=magzpt-2.5*log10(flim)                       ! 5-sigma mag limit

c *** print out photometry into ASCII catalog
      do k=1,nrows

      xcord=rec(k,3)
      ycord=rec(k,5)
      iflcol=20    ! column number of Aper_flux_1
      do iap=1,napers
      flux(iap)=max(0.1,rec(k,iflcol))
      fluxerr(iap)=max(1.0,rec(k,iflcol+1))
      iflcol=iflcol+2
      end do
      peak=rec(k,18)        ! peak flux
      peakerr=rec(k,19)     ! peak flux error
      skyloc=rec(k,56)      ! local sky
      tpeak=peak+skyloc
      confidence=rec(k,58)  ! average confidence level within default rcore aperture
      errbit=rec(k,55)      ! bit pattern listing various processing error flags
                            ! initially set to the number of bad pixels within radius 'rcore'
      ellipt=rec(k,8)
      pa=rec(k,9)-theta                                 ! correct to sky PA
      pa=90.0-pa                                        ! don't ask
c      if(pa.gt.180.0)pa=pa-180.0
      icls=rec(k,61)

c *** flag possibly saturated objects and correct flux for saturation
      ! Ring apertures are used for avoiding flux bias in PSF core:
      ! F = (f2 - f1) * C  ,
      ! C = c1c2/(c1-c2) = c1 / (c1/c2 - 1) = c1 / (delapcor - 1) ,
      ! where F is the corrected flux of a star, f1, f2 are fluxes in ap1 and ap2 (R_ap1 < R_ap2),
      ! and F = f1*c1 , F = f2*c2  (by definition)
      !
      if(tpeak.gt.saturate) then
      icls=-9                                                ! saturated ?
      ! Saturation correction using ring apertures
      ! AP7:
      fixed(7)=(flux(7)-flux(3))/(delapcor(7)-1.0)           ! flux in the ring between AP7 and AP3
      fixed(7)=fixed(7)*apcor(3)
      flux(7)=max(fixed(7),flux(7))
      ! AP6:
      fixed(6)=(flux(6)-flux(3))/(delapcor(6)-1.0)           ! flux in the ring between AP6 and AP3
      fixed(6)=fixed(6)*apcor(3)
      flux(6)=max(fixed(6),flux(6))
      ! AP5:
      fixed(5)=(flux(5)-flux(3))/(delapcor(5)-1.0)           ! flux in the ring between AP5 and AP3
      fixed(5)=fixed(5)*apcor(3)
      flux(5)=max(fixed(5),flux(3))
      ! AP4:
      fixed(4)=(flux(4)-flux(3))/(delapcor(4)-1.0)           ! flux in the ring between AP4 and AP2
      fixed(4)=fixed(4)*apcor(3)
      flux(4)=max(fixed(4),flux(4))
      ! AP3:
      fixed(3)=(flux(3)-flux(2))/(delapcor(3)-1.0)           ! flux in the ring between AP3 and AP2
      fixed(3)=fixed(3)*apcor(2)
      flux(3)=max(fixed(3),flux(3))
      ! AP2:
      fixed(2)=(flux(2)-flux(1))/(delapcor(2)-1.0)           ! flux in the ring between AP2 and AP1
      fixed(2)=fixed(2)*apcor(1)
      flux(2)=max(fixed(2),flux(2))
      ! AP1:
      fixed(1)=(flux(1)-peak)/(delapcor(1)-1.0)              ! flux in the 'ring' between AP1 and PEAK
      fixed(1)=fixed(1)*apcorpk
      flux(1)=max(fixed(1),flux(1))
      endif
c *** flag objects containing bad pixels
      if(errbit.gt.0.0)icls=-7

c *** a bit of astrometry
      xx=xcord
      yy=ycord
      call radeczp(xx,yy,xi,xn)
      call rahour(xx,ih,im,ss)
      call radeg(yy,id,in,rr)
      if(yy.lt.0)then
        csign='-'
        id=abs(id)
        in=abs(in)
        rr=abs(rr)
      else
        csign='+'
      endif

c *** apply aperture and distortion corrections
      call distort(xcord,ycord,distortcorr)                 ! field distortion correction
      do iap=1,napers
      if(icls.eq.-9) then
        flux(iap)=flux(iap)*percorr*distortcorr               ! ap. corr. already applied (see above)
      else
        flux(iap)=flux(iap)*apcor(iap)*percorr*distortcorr    ! aperture correction
      end if
      fluxerr(iap)=fluxerr(iap)*apcor(iap)*percorr*distortcorr
      magobj(iap)=magzpt-2.5*log10(flux(iap))                    ! object magnitude
      magerr(iap)=2.5*log10(1.0+fluxerr(iap)/flux(iap))
c      magerr=max(0.01,magerr)                          ! allow for systematics
      end do

c *** apply illumination correction table if required
      if(illum)then
        xi=xi*180.0/pi                                 ! xi,xn in deg
        xn=xn*180.0/pi
        do iap=1,napers
        call illumcor(tablename,xi,xn,magobj(iap))
        end do
      endif

c *** write results to output
      if(id.lt.10)then
        if(flux(3).ge.0.5)then                      ! only detected objects
        write(3,'(2i3,f8.4,1x,2a,i1,i3,f7.3,2f9.2,14f16.3,2i3,2f7.2,
     1  f10.4)')
     2  ih,im,ss,csign,czero,id,in,rr,xcord,ycord,
     3  magobj(1),magerr(1),magobj(2),magerr(2),magobj(3),magerr(3),
     4  magobj(4),magerr(4),magobj(5),magerr(5),magobj(6),magerr(6),
     5  magobj(7),magerr(7),mef,icls,ellipt,pa,confidence
        endif
      else
        if(flux(3).ge.0.5)then
        write(3,'(2i3,f8.4,1x,a,i2,i3,f7.3,2f9.2,14f16.3,2i3,2f7.2,
     1  f10.4)')
     2  ih,im,ss,csign,id,in,rr,xcord,ycord,
     3  magobj(1),magerr(1),magobj(2),magerr(2),magobj(3),magerr(3),
     4  magobj(4),magerr(4),magobj(5),magerr(5),magobj(6),magerr(6),
     5  magobj(7),magerr(7),mef,icls,ellipt,pa,confidence
        endif
      endif

      enddo                                            ! nrows loop
      enddo                                            ! mef loop

c *** and finish off
      call ftclos(iunit,status)
      if(status.ne.0)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
      endif
      end

c *** -----------------------------------------------------------------------

      subroutine fitsio_rtable(rec,iunit,mef,nrows,ncols,apcor,delapcor,
     1 apcorpk,flim,filter,exptime,utstart,dateobs,percorr,magzpt,
     2 magzrr,airmass,extinct,skyvar,gain,saturate,pltscl,theta,iesocam)
      parameter (ncmax=80,nrmax=2000000,napers=7)
      real*8 tpa,tpd,tand,secd,tpi,pi
      real*4 rec(nrmax,ncmax)
      real*4 apcor(napers)
      real*4 delapcor(napers)
      integer status,hdutype,naxes(2)
      real*4 enullval,magzpt,magzrr
      character*80 errmsg,utstart,dateobs
      character*72 comment
      character*68 value
      character*16 ttype(ncmax),tunit(ncmax),tform(ncmax)
      character*11 filter
      character*1 ciap
      logical anynull
      common /rd/ tpa,tpd,tand,secd,tpi,a,b,c,d,e,f,projp1,projp3,
     1            projp5,iwcstype
c *** FITSIO_RTABLE - reads fits binary tables
      pi=4.0*atan(1.0)
      anynull=.false.
      enullval=0.0
c *** find out what sort of table extension it is
      status=0
      call ftmahd(iunit,mef+1,hdutype,status)
      if(hdutype.eq.1)then
        print *,'Extension is an ASCII table'
        stop ' - not implemented yet'
      else if(hdutype.eq.2)then
        print *,'Extension is a BINARY table'
      endif
c *** read assorted keywords to find table dimensions etc.
      call ftgknj(iunit,'NAXIS',1,2,naxes,nfound,status)
      if(nfound.ne.2)stop 'Not enough NAXIS keywords'
      ncols=naxes(1)/4                                    ! bytes to r*4
      nrows=naxes(2)
      print *,' '
      print *,'No. of columns =',ncols
      print *,'No. of rows    =',nrows
      if(nrows.gt.nrmax)then
        print *, ' '
        print *, 'nrmax limit exceeded'
        print *, ' '
        stop
      endif
      if(iesocam.eq.0)then
        call ftgkys(iunit,'WFFBAND',filter,comment,status)
        if(status.eq.202)then
          status=0
          call ftgkys(iunit,'FILTER',filter,comment,status)
          if(status.eq.202)then
            status=0
            call ftgkys(iunit,'HIERARCH ESO INS FILT1 NAME',
     1                  filter,comment,status)
          endif
        endif
        if(status.eq.202)then
          status=0
          filter='         '
        endif
      endif

      call ftgkye(iunit,'EXPTIME',exptime,comment,status)
      if(status.eq.202)then
        status=0
        call ftgkye(iunit,'EXPOSED',exptime,comment,status)
        if(status.eq.202)then
          status=0
          call ftgkye(iunit,'EXP_TIME',exptime,comment,status)
        endif
      endif
      if(status.eq.202)then
        status=0
        exptime=0.0
      else
        exptime=abs(exptime)
        exptime=max(1.0,exptime)
      endif

      call ftgkye(iunit,'SKYLEVEL',skylevel,comment,status)
      if(status.eq.202)then
        status=0.0
        call ftgkye(iunit,'ESO DRS SKYLEVEL',skylevel,comment,status)
        if(status.eq.202)then
          status=0
          skylevel=0.0
        endif
      endif

      call ftgkye(iunit,'SKYNOISE',skynoise,comment,status)
      if(status.eq.202)then
        status=0
        call ftgkye(iunit,'ESO DRS SKYNOISE',skynoise,comment,status)
        if(status.eq.202)then
          status=0
          skynoise=1.0
        endif
      endif

      call ftgkye(iunit,'RCORE',rcore,comment,status)
      if(status.eq.202)then
        status=0
        call ftgkye(iunit,'ESO DRS RCORE',rcore,comment,status)
        if(status.eq.202)then
          status=0
          rcore=1.0
        endif
      endif

      flim=5.0*sqrt(pi*rcore**2)*skynoise           ! 5-sigma flux limit
      skyvar=(pi*rcore**2)*skynoise**2              ! for error analysis

      call ftgkye(iunit,'GAIN',gain,comment,status)
      if(status.eq.202)then
        status=0
        gain=2.0                                    ! default
      endif
      if(iesocam.eq.1)gain=4.0
      if(iesocam.eq.2)gain=2.0

      call ftgkye(iunit,'SATURATE',saturate,comment,status)
      if(status.eq.202)then
        status=0
        saturate=30000.0                            ! safe guess
      endif
      print *,'SATURATE =',saturate                 ! average saturation level in frame
      saturate=0.9*saturate

      if(status.ne.0)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
      endif

      call ftgkye(iunit,'APCORPK',apcorpk,comment,status)
      if(status.eq.202)then
          status=0
          apcorpk=0.0  ! default is no aperture correction
      endif
      apcorpk=10.0d0**(0.4d0*apcorpk)

      do iap=1,napers
        write(ciap,'(i1)') iap
        call ftgkye(iunit,'APCOR'//ciap,apcor(iap),comment,status)
        if(status.eq.202)then
          status=0
          apcor(iap)=0.0  ! default is no aperture correction
        endif
        print *,'APCOR'//ciap,'=',apcor(iap)
      end do
      ! Aperture corrections in flux for ring apertures
      delapcor(1)=10.0d0**(0.4d0*(apcorpk -apcor(1)))  ! 'ring' between ap1 and peak flux
      delapcor(2)=10.0d0**(0.4d0*(apcor(1)-apcor(2)))  ! ring between ap2 and ap1
      delapcor(3)=10.0d0**(0.4d0*(apcor(2)-apcor(3)))  ! ring between ap3 and ap2
      delapcor(4)=10.0d0**(0.4d0*(apcor(3)-apcor(4)))  ! ring between ap3 and ap4
      delapcor(5)=10.0d0**(0.4d0*(apcor(3)-apcor(5)))  ! ring between ap3 and ap5
      delapcor(6)=10.0d0**(0.4d0*(apcor(3)-apcor(6)))  ! ring between ap3 and ap6
      delapcor(7)=10.0d0**(0.4d0*(apcor(3)-apcor(7)))  ! ring between ap3 and ap7
      ! Aperture corrections in flux for apertures 1-7
      do iap=1,napers
        apcor(iap)=10.0d0**(0.4d0*apcor(iap))
      end do

      call ftgkye(iunit,'PERCORR',percorr,comment,status)
      if(status.eq.202)then
        status=0
        percorr=0.0
      endif
      percorr=10.0**(0.4*percorr)

      flim=flim*apcor(3)*percorr                    ! including ap and per corr

      if(iesocam.eq.0)then
        call ftgkye(iunit,'AIRMASS',airmass,comment,status)
        if(status.eq.202)then
          status=0
          call ftgkye(iunit,'AMSTART',airmass1,comment,status)
          call ftgkye(iunit,'AMEND',airmass2,comment,status)
          if(status.eq.202)then
            status=0
            airmass=1.0
          else
            airmass=0.5*(airmass1+airmass2)
          endif
        endif
      endif

      call ftgkye(iunit,'MAGZPT',magzpt,comment,status)
      if(status.eq.202)then
        status=0
        magzpt=-1.0
      endif

      call ftgkye(iunit,'MAGZRR',magzrr,comment,status)
      if(status.eq.202)then
        status=0
        magzrr=-1.0
      endif

c *** get extinction value - if not present use default
      call ftgkye(iunit,'EXTINCT',extinct,comment,status)
      if(status.eq.202)then
        status=0
        extinct=0.05
      endif

c *** check if WCS information present
      call ftgkey(iunit,'CTYPE1',value,comment,status)
      if(status.eq.202)then
        status=0
        call ftgkey(iunit,'TCTYP3',value,comment,status)
      endif
      if(status.ne.202)then
        if(index(value,'RA---').gt.0)then           ! need a better test
          iwcs=1
        else
          iwcs=0
        endif
      else
        iwcs=0
        status=0
      endif
      if(iwcs.eq.1)then
        print *,'Valid WCS present'
      else
        print *,'No WCS present'
      endif
c *** now get celestial coordinate information the hard way
      if(iwcs.eq.1)then
        iwcstype=1                                             ! ZPN default
        if(iesocam.eq.0)then
          call ftgkys(iunit,'CTYPE1',value,comment,status)
          if(index(value,'RA---TAN').gt.0)iwcstype=0           ! TAN only other
          call ftgkye(iunit,'CRPIX1',c,comment,status)
          call ftgkye(iunit,'CRPIX2',f,comment,status)
          call ftgkyd(iunit,'CRVAL1',tpa,comment,status)
          call ftgkyd(iunit,'CRVAL2',tpd,comment,status)
        else
          call ftgkys(iunit,'TCTYP3',value,comment,status)
          if(index(value,'RA---TAN').gt.0)iwcstype=0           ! TAN only other
          call ftgkye(iunit,'TCRPX3',c,comment,status)
          call ftgkye(iunit,'TCRPX5',f,comment,status)
          call ftgkyd(iunit,'TCRVL3',tpa,comment,status)
          call ftgkyd(iunit,'TCRVL5',tpd,comment,status)
        endif
        tpa=tpa*pi/180.0
        tpd=tpd*pi/180.0
        tand=tan(tpd)
        secd=1.0/cos(tpd)
        call ftgkye(iunit,'CD1_1',a,comment,status)
        if(status.eq.202)then
          status=0
          call ftgkye(iunit,'TC3_3',a,comment,status)
        endif
        a=a*pi/180.0
        call ftgkye(iunit,'CD2_2',e,comment,status)
        if(status.eq.202)then
          status=0
          call ftgkye(iunit,'TC5_5',e,comment,status)
        endif
        e=e*pi/180.0
        call ftgkye(iunit,'CD1_2',b,comment,status)
        if(status.eq.202)then
          status=0
          call ftgkye(iunit,'TC3_5',b,comment,status)
          if(status.eq.202)then
            status=0
            b=0.0
          endif
        endif
        b=b*pi/180.0
        call ftgkye(iunit,'CD2_1',d,comment,status)
        if(status.eq.202)then
          status=0
          call ftgkye(iunit,'TC5_3',d,comment,status)
          if(status.eq.202)then
            status=0
            d=0.0
          endif
        endif
        d=d*pi/180.0
        if(iverbose.eq.1)then
          print *,' '
          print *,'Frame transform constants:'
          print *,a,b,c
          print *,d,e,f
          print *,'Tangent points:',tpa,tpd
        endif
c *** pixel scale and orientation
        pltscl=sqrt(0.5*(a**2+b**2+d**2+e**2))*180.0/pi               ! deg/pix
        theta=0.5*(atan(abs(b)/abs(a))+atan(abs(d)/abs(e)))*180.0/pi  ! deg
c *** read zp constants from header
        projp1=1.0
        call ftgkye(iunit,'PV2_3',projp3,comment,status)
        if(status.eq.202)then
          status=0
          call ftgkye(iunit,'PROJP3',projp3,comment,status)
          if(status.eq.202)then
            status=0
            call ftgkye(iunit,'TV5_3',projp3,comment,status)
            if(status.eq.202)then
              status=0
              projp3=0.0
            endif
          endif
        endif
        call ftgkye(iunit,'PV2_5',projp5,comment,status)
        if(status.eq.202)then
          status=0
          call ftgkye(iunit,'TV5_5',projp5,comment,status)
          if(status.eq.202)then
            status=0
            projp5=0.0
          endif
        endif
      endif
c *** observation time and date
      if(iesocam.eq.0)then
        call ftgcrd(iunit,'UTSTART',utstart,status)
        if(status.eq.202)then
          status=0
          call ftgcrd(iunit,'UTC-OBS',utstart,status)
        endif
        if(status.eq.202)then
          utstart='UTST    ='
          status=0
        endif
        call ftgcrd(iunit,'DATE-OBS',dateobs,status)
        if(status.eq.202)then
          dateobs='DATE    ='
          status=0
        endif
      endif
      if(status.ne.0)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
      endif
c *** now get rest of table info
      call ftgkns(iunit,'TTYPE',1,ncmax,ttype,nfound,status)
      call ftgkns(iunit,'TUNIT',1,ncmax,tunit,nfound,status)
      call ftgkns(iunit,'TFORM',1,ncmax,tform,nfound,status)
      if(ncols.ne.nfound)stop 'Something screwy with header'
      if(iverbose.eq.1)then
        print *,' '
        print *,'Parameters in table:'
        print *,(ttype(i),i=1,ncols)
        print *,' '
      endif
c *** now read data column by column (yuk!)
      do j=1,ncols
      call ftgcve(iunit,j,1,1,nrows,enullval,rec(1,j),anynull,status)
      enddo
c *** check for errors
      if(status.ne.0)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
      endif
      return
      end

c *** -----------------------------------------------------------------------

      subroutine radeczp(x,y,xi,xn)
      real*8 tpa,tpd,tand,secd,tpi
      real*8 xi,xn,aa,alpha,delta,x,y
      common /rd/ tpa,tpd,tand,secd,tpi,a,b,c,d,e,f,projp1,projp3,
     1            projp5,iwcstype
c *** RADECZP converts x,y coordinates to ra and dec using ZP wrt optical axis
c *** NB. This is an extension of the ARC projection using Zenithal polynomials
      x=x-c
      y=y-f
      xi=a*x+b*y
      xn=d*x+e*y
      r=sqrt(xi**2+xn**2)
      if(iwcstype.eq.1)then
        rfac=projp1+projp3*r**2+projp5*r**4    ! NB this is a 1st order approx
        r=r/rfac
        rfac=projp1+projp3*r**2+projp5*r**4    ! now 2nd order correction
      else
        if(r.eq.0.0)then
          rfac=1.0
        else
          rfac=r/tan(r)
        endif
      endif
      xi=xi/rfac
      xn=xn/rfac
      aa=atan(xi*secd/(1.0-xn*tand))
      alpha=aa+tpa
      delta=atan((xn+tand)*sin(aa)/(xi*secd))
      x=alpha
      y=delta
      if(x.gt.tpi)x=x-tpi
      if(x.lt.0.0)x=x+tpi
      return
      end

c *** -----------------------------------------------------------------------

      subroutine distort(x,y,distortcorr)
      real*8 tpa,tpd,tand,secd,tpi
      common /rd/ tpa,tpd,tand,secd,tpi,a,b,c,d,e,f,projp1,projp3,
     1            projp5,iwcstype
c *** DISTORT converts x,y coordinates to standard coordinates
c ***         and works out flux distortion factor
c *** NB. This is an extension of the ARC projection using Zenithal polynomials
      xi=a*(x-c)+b*(y-f)
      xn=d*(x-c)+e*(y-f)
      r=sqrt(xi**2+xn**2)
      if(iwcstype.eq.1)then
        distortcorr=1.0+3.0*projp3*r**2/projp1+
     1                  5.0*projp5*r**4/projp1 ! this is only a 1st order app
        distortcorr=distortcorr*(1.0+(projp3*r**2+projp5*r**4)/projp1)
      else
        if(r.eq.0.0)then
          distortcorr=1.0
        else
          distortcorr=1.0/(1.0+tan(r)**2)
          distortcorr=distortcorr*atan(r)/r
        endif
      endif
      distortcorr=1.0/distortcorr              ! is this the right sign
      return                                   ! r.dr.dtheta --> r'dr'dtheta
      end

c *** -----------------------------------------------------------------------

      subroutine clobber(iunit,name)
c *** CLOBBER checks to see if file name exists and then deletes it
      integer ier,stat,statb(13),status,system
      character*(*) name
      character*80 filename
      character*3 del /'rm '/
      character*1 space /' '/
      filename=name
      i=1
      do while (filename(i:i).ne.' '.and.i.lt.80)
      i=i+1
      enddo
      open(unit=iunit,file=filename)
      ier=stat(filename,statb)
      if(statb(8).eq.0.or.ier.ne.0)then
        status=system(del//filename)
      else
        print *,'File: '//filename(1:i)//' already exists - clobbering'
        status=system(del//filename)
      endif
      close(unit=iunit)
      return
      end

c *** ------------------------------------------------------------------------

      subroutine rahour(radian,ihour,min,sec)
      real*8 radian,temp
      character*6 dum
c *** RAHOUR  converts radians into hours,mins,secs
      tpi=8.0*atan(1.0)
      if(radian.lt.0.0)radian=radian+tpi
c *** convert to secs
      temp=10800.0*abs(radian)/atan(1.0)
      ir=int(temp)
      ihour=ir/3600
      min=(ir-3600*ihour)/60
      sec=(temp-ihour*3600-min*60)

c *** check for illegal sexagesimals
      write(dum,'(f6.2)') sec
      if(index(dum,'60.').ne.0)then
        sec=0.0
        min=min+1
      endif
      if(min.eq.60)then
        min=0
        ihour=ihour+1
        if(ihour.eq.24)ihour=0
      endif

      return
      end

c *** ------------------------------------------------------------------------

      subroutine radeg(radian,ideg,min,sec)
      real*8 radian,temp
      character*5 dum
c *** RADEG  converts radians into degs,mins,secs
      iflag=1
      if(radian.lt.0.0)iflag=-1
c *** convert to secs
      temp=162000.0*abs(radian)/atan(1.0)
      ir=int(temp)
      ideg=ir/3600
      min=(ir-3600*ideg)/60
      sec=(temp-ideg*3600-min*60)*iflag

c *** check for illegal sexagesimals
      write(dum,'(f5.1)') sec
      if(index(dum,'60.').ne.0)then
        sec=0.0
        min=min+1
      endif
      if(min.eq.60)then
        min=0
        ideg=ideg+1
      endif

      min=min*iflag
      ideg=ideg*iflag
      return
      end

c *** ------------------------------------------------------------------------

      subroutine illumcor(tablename,xideg,xndeg,magobj)
      real*8 xideg,xndeg
      real*4 bin(35,35)
      real*8 magobj
      character*80 tablename
      logical first
      common /illum/ first,bin
c *** ILLUMCOR does illumination correction (currently only for WFCAM)
      if(first)then
c *** read in table
        open(unit=1,file=tablename,status='old')
        read(1,*)
        do j=1,35
        do i=1,35
        read(1,*) xi,xn,bin(i,j),numb
        enddo
        enddo
        close(unit=1)
        first=.false.
      endif
c *** use input xi,xn in deg to find appropriate closest illum corr
      i=int(xideg*50.0+18.0)
      delx=xideg*50.0+18.0-i
      j=int(xndeg*50.0+18.0)
      dely=xndeg*50.0+18.0-j
c *** bilinear interpolation
      if(magobj.ne.0.0)
     1 magobj=magobj+(1.0-dely)*((1.0-delx)*bin(i,j)+delx*bin(i+1,j))+
     2                   dely*((1.0-delx)*bin(i,j+1)+delx*bin(i+1,j+1))
      return
      end

max_stars_repo_name: carpyncho/merlomatch

max_issues_repo_name: carpyncho/merlomatch

max_forks_repo_name: carpyncho/merlomatch


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE DBPTIM (OPTION, NSTEPS, TIMES)
C=======================================================================

C   --*** DBPTIM *** (EXOLIB) Print database steps and min/max times
C   --   Written by Amy Gilkey - revised 12/18/87
C   --
C   --DBPTIM displays the number of time steps and the minimum and maximum
C   --time on the database.
C   --
C   --Parameters:
C   --   OPTION - IN - '*' to print all, else print options:
C   --                 'N' to print the number of time steps
C   --                 'M' to print the minimum and maximum time step times
C   --                 'T' to print the time step times
C   --   NSTEPS - IN - the number of time steps
C   --   TIMES  - IN - the database time steps

C   --Routines Called:
C   --   INTSTR - (STRLIB) Convert integers to strings
C   --   NUMSTR - (STRLIB) Convert numbers to engineering notation
C   --   SQZSTR - (STRLIB) Delete extra blanks from string

      CHARACTER*(*) OPTION
      INTEGER       NSTEPS
      REAL          TIMES(*)
      CHARACTER*80  STRING
      CHARACTER*5   ISTRA
      CHARACTER*20  RSTR(3)
      REAL          RNUM(3)
C     Logical flags for options
      LOGICAL ALL

      ALL = OPTION .EQ. '*'
      WRITE (*, *)

      IF (ALL .OR. (INDEX (OPTION, 'N') .GT. 0)) THEN
         WRITE (STRING, 10000) NSTEPS
10000    FORMAT ('Number of time steps = ', I8)
         CALL SQZSTR (STRING, LSTR)
         WRITE (*, 10030) STRING(:LSTR)
      END IF

      IF (ALL .OR. (INDEX (OPTION, 'M') .GT. 0)
     &    .AND. (NSTEPS .GT. 0)) THEN
         CALL MINMAX (NSTEPS, TIMES, TIMMIN, TIMMAX)

         RNUM(1) = TIMMIN
         RNUM(2) = TIMMAX
         CALL NUMSTR (2, 4, RNUM, RSTR, LSTR)
         WRITE (*, 10030)
     &      '   Minimum time = ', RSTR(1)(:LSTR)
         WRITE (*, 10030)
     &      '   Maximum time = ', RSTR(2)(:LSTR)
      END IF

      IF (ALL .OR. (INDEX (OPTION, 'T') .GT. 0)
     &    .AND. (NSTEPS .GT. 0)) THEN
         CALL INTSTR (1, 0, NSTEPS, ISTRA, LISTR)
         RNUM(2) = TIMES(1)
         IF ((RNUM(2) .EQ. 0.0) .AND. (NSTEPS .GT. 1))
     &      RNUM(2) = TIMES(2)
         RNUM(3) = TIMES(NSTEPS)
         DO 110 I = 1, NSTEPS
            CALL INTSTR (1, LISTR, I, ISTRA, LI)
            RNUM(1) = TIMES(I)
            CALL NUMSTR (3, 4, RNUM, RSTR, LSTR)
            WRITE (*, 10020) ISTRA(:LI), RSTR(1)(:LSTR)
10020       FORMAT (1X, 3X, 'Step ', A, ')', 3X, A, :, 3X, A)
  110       CONTINUE
      END IF

      RETURN
10030  FORMAT (1X, 5A)
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute mean 2PCF multipoles from Nres=1000 EZmocks
! for any issue please contact ginevra.favole@port.ac.uk
!-------------------------------------------------------
!compile as:
!  ifort -fopenmp mean_multipoles_rsd.f90 -o meanmrsd
!execute as:
! ./meanmrsd

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program main
implicit none
CHARACTER(len=150) :: FN
CHARACTER(len=150) :: GN
CHARACTER(150), PARAMETER   :: output='../multipoles_EZmocks3Gpc/multipoles_5Mpc/mean_multipoles_rsd.txt'

integer, parameter :: Nres=999,cov_num=40
integer :: i, j, k, io, N_halsum
real*8 :: mxi0(1:cov_num), mxi2(1:cov_num), mxi4(1:cov_num), s(1:cov_num), temp
real*8, dimension(:), allocatable  :: bin_mock, xi0_mock, xi2_mock, xi4_mock

!------------------------------------------------ format input read all
mxi0=0.0
mxi2=0.0
mxi4=0.0

DO I=0,Nres
WRITE(FN,10)I
WRITE(6,*)FN
j=1
OPEN(1,FILE=FN)
DO
READ(1,*,IOSTAT=io) temp
IF (io < 0) THEN
N_halsum = j-1
EXIT
ENDIF
j =j+1
ENDDO

allocate(bin_mock(1:N_halsum),xi0_mock(1:N_halsum),xi2_mock(1:N_halsum),xi4_mock(1:N_halsum))
rewind(1)
do k=1,N_halsum
read(1,*,IOSTAT=io) bin_mock(k),xi0_mock(k),xi2_mock(k),xi4_mock(k)
enddo

do k=1,N_halsum
s(k) = bin_mock(k)
mxi0(k) = mxi0(k)+xi0_mock(k)
mxi2(k) = mxi2(k)+xi2_mock(k)
mxi4(k) = mxi4(k)+xi4_mock(k)
enddo
CLOSE(1)
DEALLOCATE(bin_mock, xi0_mock, xi2_mock, xi4_mock)
ENDDO

open(unit=1020,file=output)
do k=1,N_halsum
mxi0(k) = mxi0(k)/dble(Nres+1.)
mxi2(k) = mxi2(k)/dble(Nres+1.)
mxi4(k) = mxi4(k)/dble(Nres+1.)
write(1020,'(f10.4,x,f20.8,x,f20.8,x,f20.8)') s(k), mxi0(k), mxi2(k), mxi4(k)
enddo
close(1020)

10 FORMAT('../multipoles_EZmocks3Gpc/multipoles_5Mpc/xilrsd_',I3.3,'.txt')
end program